Mega‐Tidal and Surface Flooding Controls on Coastal Groundwater and Saltwater Intrusion Within Agricultural Dikelands

Climate change will increase sea levels, driving saltwater into coastal aquifers and impacting coastal communities and land use viability. Coastal aquifers are also impacted by tides that control groundwater‐ocean interactions and maintain an “upper saline plume” (USP) of brackish groundwater. Coastal dikes are designed to limit the surface impacts of high‐amplitude tides, but, due to ongoing sea‐level rise (SLR), low‐lying dikelands and underlying aquifers are becoming increasingly vulnerable to flooding from high tides and storm surges. This study combines field observations with numerical modeling to investigate ocean‐aquifer mixing and future saltwater intrusion dynamics in a mega‐tidal (tidal range >8 m) dikeland along the Bay of Fundy in Atlantic Canada. Field data revealed strong connectivity between the ocean and coastal aquifer, as evidenced by pronounced tidal oscillations in deeper groundwater heads and an order of magnitude intra‐tidal change in subsurface electrical resistivity. Numerical model results indicate that SLR and surges will force the migration of the USP landward, amplifying salinization of freshwater resources. Simulated storm surges can overtop the dike, contaminating agricultural soils. The presence of dikes decreased salinization under low surge scenarios, but increased salinization under larger overtopping scenarios due to landward ponding of seawater behind the dike. Mega‐tidal conditions maintain a large USP and impact aquifer freshening rates. Results highlight the vulnerability of terrestrial soil landscapes and freshwater resources to climate change and suggest that the subsurface impacts of dike management decisions should be considered in addition to protection measures associated with surface saltwater intrusion processes.

Coastal aquifer dynamics and the width of the transition zone between salt-and freshwater are strongly influenced by inland and tidal forcing (Heiss & Michael, 2014;LeRoux et al., 2021;Pool et al., 2014;Robinson et al., 2007;Robinson & Li, 2004).Analytical solutions have been developed to quantify how tidal signals propagate into coastal aquifers and influence groundwater dynamics (Bakker, 2019;Inouchi et al., 1990;Li & Jiao, 2001;Nielsen, 1990).Such solutions have been applied in an inverse manner to infer aquifer properties from observed groundwater tidal signals that are damped and lagged in comparison to surface water tides (Zhang et al., 2020).Tidal pumping and density differences also drive the formation of the upper saline plume (USP, Figure 1a), which is a brackish zone perched above the salt wedge that is a "hot spot" for biogeochemical reactions (Robinson  et al., 2018).The presence and morphology of USPs are strongly dependent on factors including beach morphology, tidal range, density differences, and heterogeneity (Bakhtyar et al., 2013;Evans & Wilson, 2016;Robinson et al., 2018).The USP dimensions and spatial and temporal stability are also dependent on the magnitude and timing of external forcing (Buquet et al., 2016;Fang et al., 2021;Fang, Zheng, Guo, et al., 2022;Fang, Zheng, Wang et al., 2022;Grünenbaum et al., 2023;Kuan et al., 2012;Yu et al., 2019).Despite the relationship between tidal amplitude and USP morphology and dynamics, numerical modeling studies investigating USP dynamics, or more generally aquifer-ocean interactions, have been limited to micro-(and occasionally macro-) tidal ranges (e.g., Levanon et al., 2017;Yang et al., 2013).To our knowledge, there are no published studies of field and modeling investigations of ocean-aquifer mixing in mega-tidal settings, where tidal ranges exceed 8 m.
Dikes along coastlines are built to disrupt the natural surface water tidal regime and provide the first line of defense during high tides or coastal storms (Rahman et al., 2019;Reed, 1995;Richards & Daigle, 2011).While hard sea barriers like dikes can prevent saltwater flooding, they also result in the loss of salt marshes that dissipate wave energy and provide important wildlife habitat (Möller & Spencer, 2002;Narayan et al., 2017;Sherren et al., 2021;White & Kaplan, 2017).Coastal dikes in many regions of the world were not initially designed for SLR or the greater intensity storms projected to occur under a changing climate.Consequently, many of these barriers are at risk of being breached, allowing saltwater to flood land behind dikes and salinize freshwater resources and terrestrial environments such as crop land (Figure 1d) (Carrión-Mero et al., 2021;Elsayed & Oumeraci, 2018a).In addition to loss of suitable crop land, there are ecological consequences of SWI including biogeochemical changes in coastal substrates, loss of coastal forests, freshwater wetland conversion, expansion of invasive species, and declines in agricultural productivity (Briggs et al., 2021;Brinson et al., 1995;Guimond & Michael, 2020;Tully et al., 2019).Many knowledge gaps exist related to the future fate of dikelands and the ecosystem services they provide (Sherren et al., 2021).Potential dike management solutions to the present and intensifying threat of climate change include increasing dike elevations, retreating dikes back to increase the foreshore marsh width, or removing dikes altogether (van Proosdij et al., 2010).While the surface impacts of these dike management solutions have been investigated, little to no research has been conducted on the impacts of dike alterations on coastal groundwater resources.
The purpose of this study is to investigate ocean-aquifer exchanges in settings with tidal ranges larger than previously considered, and to assess how coastal aquifers are influenced by dikes under present and future climate conditions.We collected hydrogeological and geophysical field data and conducted integrated 2D physics-based coupled surface and subsurface modeling to specifically address four key objectives: (a) investigate how mega tides drive coastal aquifer mixing and groundwater flow dynamics, (b) simulate how SLR will influence these dynamics, with a focus on USP migration, (c) simulate how future overtopping events during storm surges could salinize agricultural land and potentially freshwater aquifers, and (d) determine how these dynamics would change with and without dikes or tides.Ultimately, our goal is to provide new hydrogeological understanding of coastal groundwater dynamics in dikelands to help inform optimal dike management alternatives in the face of climate change.

Study Site
The Bay of Fundy in southeastern Canada and the northeastern United States (Figure 2) has some of the largest tidal ranges in the world (Garrett, 1972).Dikes built in the 1700s were developed to protect burgeoning agricultural land from tidal influences (Desplanque & Mossman, 2004;Sherren et al., 2021).The Annapolis Valley (Figure 2b) in Nova Scotia, Canada lies along the Bay of Fundy and contains many towns and field croplands in close proximity to the Fundy coast (Rivard et al., 2006).The town of Wolfville (Figure 2c, population of ∼4,500) contains several kilometers of dikes that protect terrestrial land from SWI (Sherren et al., 2021).The town and surrounding region are particularly vulnerable to coastal changes and resultant aquifer responses because they have critical infrastructure and agricultural fields (e.g., vineyards) within 100 m of the coastline (Figure 2c), and certain dikes in the area are already experiencing overtopping during spring high tides (van Proosdij & Page, 2015).Also, the town relies on fresh groundwater for domestic and industrial use, and several coastal wells have already been abandoned due to salinity.
The climate in the Annapolis Valley is relatively humid and temperate, with 1,181 mm of mean annual  precipitation and a mean air temperature of −4°C in the winter months (January to March) and 18°C in the summer months (June to August) (Environment and Climate Change Canada, 2023).SLR along the Fundy coastline is projected to be among the highest in Canada (e.g., up to 1.5 m by 2100, James et al., 2021), which increases the vulnerability of the dikes to climate change when considered in conjunction with the potential for concurrent storm surges and spring tides (e.g., McLaughlin et al., 2022).Storms in the Bay of Fundy are primarily driven by Atlantic winter storms or hurricanes that are less impactful by the time they reach higher latitudes like Atlantic Canada (e.g., McLaughlin et al., 2022).Storm surges present the greatest flooding risks when they are coincident with high tides, as typical significant wave heights and surges are on the order of 1.5 m during storms (Mulligan et al., 2019) and thus much lower than the tidal amplitude.Historically, storm tides (coincident high tides and surges) have wreaked havoc in this region, with the most devastating storm tide, the 1869 Saxby Gale, causing severe flooding and damage throughout the Fundy region (Desplanque & Mossman, 1999).The tides in the Bay of Fundy are dominated by the principal lunar tidal constituent (M2), with the tidal range varying between 5 and 16 m (moving southwest to northeast) due to the natural resonance of the bay.Given the high tidal ranges, a storm surge as low as 1.05 m, if coincident with spring high tide, would overtop the maximum dike heights (McLaughlin et al., 2022).
The bedrock geology in the Fundy region largely determines the availability of groundwater in deeper aquifers and may influence patterns of relic salinity (e.g., Seibert et al., 2023) from tidal flooding in the pre-dike era.Surficial geology affects shallow groundwater resources and soil processes.The surficial geology is characterized 10.1029/2023WR035054 5 of 21 by glacial lake deposits (sand, silt, clay), glacial river deposits (sand and gravel), glacial till, and marine deposits (Rivard et al., 2012).The bedrock geology is characterized by thickly bedded, fine-to coarse-grained sandstone (Wolfville formation), and is considered one of the most important aquifers in the region (Rivard et al., 2012).Groundwater can be extracted in relatively high quantities (e.g., average yield of 34 L/s per well, Rivard et al., 2012) in the study area, particularly in the Wolfville formation or glaciofluvial sand and gravel deposits, both of which can be found where most of the population is located (Rivard et al., 2012).The town of Wolfville uses the aquifer in these surficial sediments as part of their water supply (Hennigar & Kennedy, 2006;Rivard et al., 2006).Understanding the groundwater-basin interactions and saltwater flushing dynamics is critical as fresh groundwater is used extensively for industrial and domestic drinking water supply (Sherren et al., 2021), and there is widespread concern about the impacts of climate change on the dikelands in this region, which is now a designated UNESCO World Heritage Site (George, 2013).

Field Instrumentation
Field methods were used to characterize the hydrodynamic and hydrogeologic conditions of the field site and to help constrain a numerical model as described in Section 3.2.To understand present-day marine drivers of groundwater dynamics, a tidal station and wave buoy were installed, a deep drilled well was instrumented, and geophysical data were collected (Figure 2c).As no government tidal station is remotely close to the study site, a tidal station (Onset HOBO, Massachusetts, level logger model U20-001-01-Ti, ±1 cm accuracy recording at 15-min intervals) was affixed to rebar and a large weight and deployed in the Minas Basin from September 2020 to June 2021 to capture the surface water tidal dynamics.These total pressure data were corrected with atmospheric pressure from a nearby barometric logger (Figure 2c, Solinst, Ontario, Barologger model, ±0.05 kPa) to yield tidal water levels.A wave buoy (SOFAR Spotter, San Francisco, frequency of 2.5 Hz recording at 30-min intervals, ±2 cm accuracy, Raghukumar et al., 2019) was installed in the basin (Figure 2c) to characterize the wave dynamics from November 2020 to November 2021 and determine if waves should be considered in the numerical model (Section 3.2).
An abandoned monitoring well near the sewage treatment plant (Figure 2c) was instrumented from November 2020 to August 2021 with a level-temperature-conductivity (LTC, Solinst Ontario, Levelogger 5 LTC, ±1 cm) logger to monitor groundwater levels and electrical conductivity (EC) in the deep bedrock aquifer, and investigate aquifer tidal oscillations.This well is located 100 m inland from the mean tide mark and is screened at a depth of 70 m below the ground surface (Figure 2c).Frequency domain analyses for the tidal and groundwater levels were conducted in Matlab using the Fast Fourier Tranform function (Frigo & Johnson, 1998) to obtain the dominant tidal constituents.Time-lapse, time-domain electromagnetic geophysical measurements of the vertical bulk resistivity profile were recorded approximately hourly in the same location in front of the dikes (seaward side, Figure 2c) from low tide to spring-high tide on 25 June 2021.Five soundings were collected using the ABEM WalkTEM (Guideline Geo).A 20 × 20 m transmitter loop was set up with short pulses of electricity, applying a 10 ms dual-moment (8A high-moment, 1A low-moment), with a powerline frequency of 50 Hz, to characterize the subsurface resistivity (e.g., Alexopoulos et al., 2017;Briggs et al., 2021;Kirsch, 2006).Aarhus University HydroGeophysics Group software "SPIA" was used to invert soundings to create resistivity models with a standard depth of investigation and yield the bulk resistivity of the matrix.This survey was used to gain a first-order understanding of the USP configuration and assess the magnitude of the intra-tidal resistivity changes to yield insight into subsurface tidal pumping.

Numerical Simulations
The numerical modeling was not conducted in an attempt to exactly reproduce aquifer conditions at the site but rather to investigate aquifer-basin interactions in mega-tidal settings and the impacts of oceanic climate change on these dynamics.More specifically, numerical model simulations were performed to investigate the present-day coastal aquifer and USP dynamics and the impacts of projected climate change on subsurface salt distributions.The water flow and solute transport simulations were performed in HydroGeoSphere, a three-dimensional, fully integrated (coupled surface-subsurface), control volume finite element model that utilizes parallel processing and adaptive time-stepping (Aquanty, 2022;Hwang et al., 2014;Therrien et al., 2010).This model was selected for its ability to simulate surface flow dynamics during storm surges, along with variable-density, variable-saturation subsurface flow.HydroGeoSphere is increasingly used for coastal groundwater applications to consider perturbations from SLR (Yang et al., 2015a), surges (Guimond & Michael, 2020;Paldor & Michael, 2021;Yu et al., 2016) and tides (Yang et al., 2015b).HydroGeoSphere solves the two-dimensional diffusion wave equation for surface flow and the Richards equation for three-dimensional variably saturated subsurface flow (Aquanty, 2022;Therrien et al., 2010).Variable density influence on subsurface flow is represented with the equivalent freshwater head approach (Frind, 1982) that uses a dimensionless relative density based on salt concentration.Surface and subsurface solute (salt) transport was simulated with consideration of both advection and dispersion.Full governing equations for flow and transport along with the surface-subsurface coupling approach are presented in Aquanty ( 2022).
The numerical model domain is 1.4 km long and 80 m deep and represents a two-dimensional profile extending between the Minas Basin and the town of Wolfville (Figures 2c and 3a).Model topography was defined using a simplified version of the provincial digital elevation model (Natural Resources and Renewables, 2006) (Figure 3a).The subsurface model domain consists of 3D quasi-rectangular, unit-thickness elements that are 1.0-2.5 m in height, and 2.5-10 m in width.Smaller elements are located within the dynamic intertidal region (Figure 3a, x = 550-1,050 m), and the larger elements are in the landward region where the tidal influence is lower.While the mesh in this study has coarser discretization than meshes typically used for simulations related to vadose zone analysis, the focus herein was on the impacts of SLR, overtopping, and tidal forcing on SWI that occurs throughout the aquifer.The total number of elements and nodes were 10,955 and 22,608, respectively.Element size and mesh discretization are in line with other HydroGeoSphere-based studies that involve simulation of variable-density salt transport and unsaturated zone dynamics (Guimond & Michael, 2020;Paldor & Michael, 2021;Yang et al., 2015a;Yu et al., 2016).In relation to the aforementioned studies, the simulations that underpin this study have the added complexity of mega-tidal boundary conditions.Due to the number of elements, the non-linear coupled equations, the short time steps, and the tight convergence criteria, the simulation wall-clock times extended up to several weeks for a 1-year simulation that utilized multiple cores of a Xeon Gold Boundary conditions were defined for both the surface and subsurface domain boundaries (Figure 3a).The surface domain extends across the top of the subsurface domain (Figure 3a) where, on the terrestrial region (x = 0-700 m, Figure 3a), a specified freshwater recharge of 100 mm yr −1 was applied as a representative value based on the previous Wolfville groundwater modeling study of Hennigar and Kennedy (2006).This recharge is also generally in line, albeit a bit lower, with baseflow-derived estimates of recharge for the broader watershed (Figure 4 of Kennedy et al., 2010).On the subtidal region of the surface domain (x = 700-1,400 m, Figure 3a), a time-varying head with a range of 12 m and seawater concentration (C/Co = 1) were specified.This head range represents an average of the tidal range observed at the tidal monitoring station as described later.A simple harmonic function was used to yield a repeating sinusoidal boundary condition that resolved the tidal cycle in five-minute increments during peaks and troughs, and up to 6.25-hr increments during the flood and ebb tides.Although incorporating spring and neap tidal cycles has been proven to alter the USP configuration (Buquet et al., 2016), incorporating these lower frequency tidal constituents did not substantially impact our modeled USP as described in Section 4.3.1.The vertical landward boundary condition was a specified freshwater head of 1 m above sea level (C/Co = 0), with the boundary location and value based on previous modeling investigations in the region and the availability of other well data (Figure 3a) (Hennigar & Kennedy, 2006).The bottom and vertical seaside boundaries were left  3b.

Table 1 Model Simulation Identifiers (ID's), Simulation Types, and Detailed Descriptions
as default, no-flow and no-solute transport boundaries (Figure 3a).The location of the seaside vertical face was based on the estimated offshore groundwater divide location inferred from the low point (i.e., in the inflowing Cornwallis River channel, Figure 2c) in the local bathymetric contour data (Figure 3a).Different basin boundary locations were tested and were found to have little impact on the flow dynamics.The bottom boundary was positioned deep into the bedrock to eliminate surface boundary effects in the shallow aquifer.
As summarized in Figures 3b and Table 1, a sequence of model runs was undertaken, and complexity was gradually increased to accommodate the complex forcing and dynamics associated with the simulations.Steady-state simulations were run using a constant mean sea-level boundary condition to achieve an equilibrium flow and transport solution that could be used as initial conditions for transient simulations (Figure 3b, run 1, Table 1).
Upon achieving an equilibrium, this simulation was used as an initial condition for transient simulations with tides turned on, for which parameters (e.g., hydraulic conductivity, specific storage) were adjusted to yield results that provided simulated heads in agreement with observation well data (Figure 3b, run 2, Table 1).Once a dynamic equilibrium was achieved, variable-density effects were activated to develop the USP (run 3, Table 1).This approach limited instabilities due to mega-tidal conditions.
The variable-density simulation (run 3, Table 1) then served as the initial conditions for simulations with climate change projections to assess lateral and vertical SWI impacts (Figure 3b).SLR simulations were run using a mean sea-level boundary that gradually increased between 2022 and 2050 at a rate in accordance with the SLR projections for the representative concentration pathway 2.6 (+0.4 m or +1.43 cm/yr, Figure 3b, run 4, Table 1) based on multi-model results compiled by James et al. (2021).After the final sea level was reached, the tides were turned back on, with the "new," final 2050 sea-level incorporated (run 5, Table 1).Due to the extensive simulation times and the need to isolate the impacts of SLR versus surges, the SLR simulations were not used as initial conditions for the surge scenarios.Rather, new sea-levels for 2050 (+0.4 m) and 2100 (+0.8 m) were applied to the tidal levels following the storm events.Storm surge simulations were run using a high-tide boundary with a surge superimposed for 2 hr (+1 m surge for 2050, +2 m surge for 2100, Figure 3b, runs 6 and 7, Table 1).These surge values were based on historic storm surge events in the Bay of Fundy recorded by tide gauges and modeled storms from historic wind records (Mulligan et al., 2019;Webster et al., 2012).This duration was chosen as the approximate time that a surge could be coincident with high tide as surges at mid or low tides do not pose any overtopping risks at this site.Additionally, a simulation was run with the dike removed and the +2 m surge imposed to assess how dike removal would impact ocean surge inundation (run 8, Table 1).Following the surge events, the tide boundary returned to the regular mega-tidal oscillation, with the 2050 sea levels incorporated (runs 9, 10, 11, Table 1).Lastly, following the +2 m surge with and without the dike, these simulations were run with and without the tides (i.e., steady-state mean sea level) to assess how the mega tides influence post-surge aquifer flushing (runs 12 and 13, Table 1).
Surface and subsurface properties were based on previous studies, field investigations, default model values, and numerical modeling iterations (Table 2).The coupling between the two domains was carried out using the dual-node approach (Therrien et al., 2010), with a coupling length of 0.1 m.This value is somewhat higher than values used in other HydroGeoSphere studies (e.g., Paldor & Michael, 2021;Yang et al., 2013) to decrease simulation run times, as the mega tidal constraints imposed in this study had especially high computational rigor.
The subsurface domain was divided into three stratigraphic units-clay, sand, and bedrock (Figure 3a, Table 2), which were selected based on previous geological mapping (Rivard et al., 2012) and numerical modeling studies (Hennigar & Kennedy, 2006) in the region.The hydraulic conductivity and storage values for these units were selected (within calibration limits) based on an iterative process that was used to reproduce the tidal groundwater signal observed in the monitoring well (Section 4).Default van Genuchten (1980) parameters for HydroGeo-Sphere were used for the variably saturated clay layer (Table 2), as a detailed investigation of the unsaturated zone dynamics was beyond the scope of this study.The dispersivity values used in this study (Table 2) are representative of sandstone (Gelhar et al., 1992) and are in line with values used in other studies (e.g., Delsman et al., 2014;Oude Essink et al., 2010;Seibert et al., 2023).

Field Results
Surface water levels in the basin are characterized by a strong semi-diurnal tidal constituent, with approximately two peaks and two troughs every day (Figure 4a, light blue).Even though the low-low tides were not always recorded as the tidal station was occasionally not submerged, the peak tidal range is approximately 14 m.The frequency domain analysis shows the dominant constituent in the basin is the principal lunar M2 constituent with a period of 12.42 hr and an amplitude of 5.41 m (Figure 4b).Other tidal constituents identified in these series include the lunar fortnightly (Mf, period of 327.86 hr), the principal lunar (O1, period of 25.82 hr), the luni-solar (K1, period of 23.92 hr), the larger lunar elliptic (N2, period of 12.66 hr), and the principal solar (S2, period of 12.00 hr) (Figure 4b) (Kamphuis, 2010).The tidal form factor (Dean & Dalrymple, 2004) obtained from the data record at our tidal station was 0.05, indicating semi-diurnal dominance.The tidal ranges during spring tides can exceed the tidal range during neap tides by 5 m (Figure 4b).
The tidal forcing from the basin propagates into the bedrock aquifer, resulting in pronounced tidal groundwater signals (up to 1.5 m range or ∼10% of the surface water range) at ∼100 m inland from the coast and 70 m below the land surface based on screen depth (Figure 4a, dark blue).The frequency domain analysis yields an amplitude of 0.58 m in the aquifer (well) for the principal lunar M2 constituent (Figure 4b).This efficient tidal signal transfer into the aquifer is indicative of confined or semi-confined aquifer conditions with high hydraulic diffusivity arising from storage effects that are limited to elastic, rather than dewatering, dynamics (Jiao & Post, 2019).In contrast with the surface water and groundwater level data, the EC in the basin and aquifer do not fluctuate and are approximately 44 mS cm −1 (seawater) and 1.2 mS cm −1 (near fresh but not potable), respectively (Figure 4a).The constant aquifer salinity at the well reveals that limited advective disturbances penetrate to the well and that the level oscillations are mainly a pressure celerity effect.
Data collected from the wave buoy in the basin (Figure 2c) are shown in Supporting Information S1 (Figure S1).The average significant wave height (i.e., average amplitude of the highest 1/3 of the waves, McCoy & Bascom, 2021) is 0.13 m, with an average peak wave period of 2.6 s (after filtering).Thus, most waves measured within the basin are smaller, irregular waves, classified as ordinary gravity waves driven by wind forces (Munk, 1950).Data from the wave buoy and tidal station indicated, as expected, that the basin is tide-dominated.Thus, tides and storm surges (rather than waves) control overtopping dynamics for a given mean sea level (McLaughlin et al., 2022).While we acknowledge these field data were only collected for a short period that did not include any extreme coastal storms, the wave data results generally align with past hydrodynamic numerical modeling results showing significant wave  heights at the study site on the order of only 0.3 m even during coastal storms (Figure 7 of Mulligan et al., 2019).
The limited wave height is due to the fetch-limited conditions and wave transformation processes arising from the complex bathymetry.Accordingly, waves were not included in the marine boundary conditions in the model.Rather, the overtopping simulations relied on tidal data and SLR and storm surge projections as described in Section 3.2.
The time-lapse, single-location electrical resistivity (inverse of conductivity) data from the seaside of the dike (Figure 2c) revealed tidal dynamics in the USP (Figures 4c and 4d).The first measurement (dark blue) yielded the highest resistivity profile, as this point was taken at low tide when groundwater flow is primarily seaward, and when the tidal intrusion of high-salinity (low-resistivity) porewater is at its minimum.As the tide rises, the deeper resistivity decreases with each timestep as more saline water intrudes into the subsurface (Figure 4c).However, in the upper 20 m, the resistivity remains low and constant over the tidal cycle, indicating that the subsurface tidal pumping is less pronounced in this upper aquifer.Collectively, these data indicate pronounced aquifer-basin exchange and the presence of a large, dynamic USP.

Discussion of Field Results
The field data yielded useful insights into the local basin and aquifer properties, ocean and aquifer interactions, and processes influencing SWI.The geophysics survey highlighted how highly dynamic the deep intertidal region is during spring high tides (Figure 4c), how these dynamics create the brackish USP, and how the deeper confined aquifer conditions "push" the salt wedge offshore as observed elsewhere along Nova Scotia's coastline (e.g., Threndyle et al., 2022).Several numerical modeling studies have simulated USP dynamics (Robinson et al., 2018;Xiao, Li, et al., 2019;Xiao, Wang, et al., 2019;Yu et al., 2019), but field observations are less common (Geng & Boufadel, 2017).This is the first reported investigation of USP dynamics in mega-tidal settings which is important to investigate in these settings given the control of tides on USP morphology (Buquet et al., 2016;Fang, Zheng, Guo, et al., 2022).Given the confined conditions in the lower aquifer and resultant high hydraulic diffusivity, the tidal pressure signals can penetrate far inland (Figure 4a) and the intra-tidal resistivity variations are pronounced (Figure 4c) in the deeper subsurface (Jacob, 1950;Jiao & Post, 2019).Such dynamics exert key controls on coastal groundwater elevations (Pool et al., 2014) and salinity, hydraulic gradients (Ataie-Ashtiani et al., 1999), and resultant flow dynamics (Fang, Zheng, Guo, et al., 2022;Fang, Zheng, Wang, et al., 2022).The attenuated signal transfer in the upper aquifer or aquitard inferred from the intra-tidal stability in the electrical resistivity (Figure 4c) is indicative of low hydraulic diffusivity arising from a combination of potentially lower permeability and higher storage due to unconfined conditions.Low-permeability, unconfined aquifers are common along the hinterlands and coastline of Nova Scotia (e.g., Craddock et al., 2022) due to the almost ubiquitous layer of glacial till or tidal mud.
The high tidal range in the groundwater level (Figure 4a) at 100 m inland from the coast and 70 m depth below the ground surface, was somewhat unexpected as other studies have shown signals that are more damped and lagged at shorter distances inland from the coast (Inouchi et al., 1990;Li & Jiao, 2001;Slooten et al., 2010).In addition to the confined conditions, the efficient tidal signal transfer is indicative of high permeability from highly fractured bedrock as has been reported elsewhere for more inland geological investigations near the study site (Kennedy & Fisher, 2022).While other studies have reported tidally controlled groundwater dynamics (Fang, Zheng, Guo, et al., 2022;Levanon et al., 2017;Robinson & Li, 2004), mega-tidal forcing on groundwater systems have been less well-documented except for sediment porewater dynamics (LeRoux et al., 2021;Stark & Hay, 2014).

Model Calibration
Comparison to field data during calibration was conducted to ensure the conceptual model and parameterization (Table 2) were reasonably representative of conditions at the site.The key calibration points included the monitoring well water levels and EC, and the time-lapse geophysical survey (Figure 2c).In particular, the low EC in the monitoring well (Figure 4a) imposed a limit on the landward extent of the saltwater wedge at that depth (Figure 5b, green diamond).The saltwater-freshwater interface location (or absence) in coastal aquifers is strongly influenced by hydraulic conductivity and is a key calibration target with a high signal-to-noise ratio compared to hydraulic heads (Pavlovskii et al., 2022).The simulated and observed dynamic heads (Figure 5c) at the location of the deep well (Figure 2c) were also used to calibrate the hydraulic conductivity and storage of the geologic units in the model (Table 2).The geophysics survey suggested the presence of a brackish USP (Figure 4d) but did not identify the presence of a clear salt wedge, in agreement with the monitoring well EC data (Figure 4a).The calibrated horizontal hydraulic conductivity for clay (1E−5 m/s) is within the range of previously tabulated values for clay (Freeze & Cherry, 1979, p. 26).The calibrated horizontal bedrock hydraulic conductivity (1E−4 m/s) is high for bedrock units, but this is a highly fractured sandstone unit (dual-domain flow) that is highly conductive as evidenced by the dynamic tidal groundwater signals (Figures 4a and 4b).Also, the calibrated conductivity is in the range of hydraulic conductivity reported for fractured bedrock in Nova Scotia (Kennedy, 2022).
The simulated head oscillations at this monitoring well location generally matched the observed heads collected in the field (Figure 5c).For example, the mean tidal amplitude from the model was 2.1 masl, while the mean tidal amplitude in the observed data was 2.4 masl.However, the fit between the modeled and observed series is not exact, particularly during spring tides, as the tidal boundary condition in the model was a simplified (single harmonic) version of the true tidal series that had multiple tidal constituents and associated frequencies (Figure 4b).An additional test simulation was run with the true tidal series (i.e., spring and neap signals included), but these additional tidal constituents did not substantially change the salinity distribution (Figure S2 in Supporting Information S1), and thus were not included in all scenarios.Lastly, the cross-shore location of the USP identified in the time-lapse geophysical data was generally reproduced in the model (Figure 5b).The results from this calibration were used for the initial conditions for the SLR and storm surge simulations (run 3, Table 1).

Sea-Level Rise Impacts on the Upper Saline Plume
The salt distributions for SLR simulations 4 and 5 (Table 1) are presented in Figure 6.After only 1 year of SLR (2022-2023), the tide-and density-driven USP nearly disappeared as the tides were "turned off" during the SLR period for numerical/computational reasons as previously described (Figure 6a).After 28 years of SLR without tides (i.e., 2022-2050), the USP has completely disappeared (Figure 6b).One-year post-SLR (2050-2051), with the tides re-introduced to investigate the extent of tidal pumping and the USP for the new sea level conditions (+0.4 m), the USP has re-developed but extends deeper (up to 45 mbsl) and further landward (Figure 6c).The USP volume after SLR increased by 19% compared to run 3 after 1 year of simulations with tides introduced.These results show that SLR drives an increase in the vertical and horizontal extent of the USP, and the USP configuration is highly sensitive to the maximum sea level elevation (Figure 6c).The USP also begins to develop salinity patterns that indicate a lack of USP stability.

Surge Impacts on the Shallow Agricultural Soils and Upper Saline Plume
The salt distributions for surge simulations (runs 9 to 13, Table 1) are presented in Figure 7. Immediately following the 1-m surge, 2-m surge, and 2-m surge with the dike removed (Figures 7a, 7c, and 7e), there is practically no difference between the initial and new 1% isocontour lines that define the USP.However, the USP begins to change within months following these surge events (Figures 7b, 7d, and 7f), shifting from stable to unstable, in line with previous modeling work (Fang, Zheng, Guo, et al., 2022).USP expansion following the surge events is not due to the surge itself but rather relates to the SLR boundary condition imposed during flushing.The 1-m surge did not induce dike overtopping but increased the vertical (downward) extent of the USP by disturbing its dynamic equilibrium (Figure 7b), and the overall simulation increased the USP volume by 31%.The 2-m surge, with the dike in place, caused overtopping which resulted in shallow soil salinization, USP expansion (+37%) (Figures 7d and 7g), and an increase in the total saltwater volume by 10%, where saltwater is defined as C/Co > 0.5.For context and reference, the 2-m surge with the dike removed immediately flooded the agricultural land and resulted in USP migration and expansion (+43%) and shallow soil and groundwater salinization (Figures 7f and 7k).For this run, the saltwater volume increased by 22.1% after 1 year and 40.3% after 4 years post-surge, respectively.
By focusing on the shallow agricultural land following the 2-m surges (Figures 7g, 7h, 7k, and 7l), the results show that the infiltration for the no-dike scenario (Figures 7k and 7l) results in higher shallow salinity concentrations compared to the same surge event with the dike present (Figures 7g and 7h), as the dike delays inundation.Additionally, the dike retains water on the landward side following the surge, resulting in "pooling" on the landside of the dike that facilitates the development of a plume immediately landward of the dike (Figure 7g) due to the greater vertical hydraulic gradients.In contrast, when the dike is removed, the infiltration is more uniform across the agricultural land, with no focused plume developing behind the dike (Figure 7k).After 4 years post-surge, salt fingers (Post & Houben, 2017) developed at >15 m below the ground surface for the dike scenario (Figure 7h), while limited salt fingers occurred for the no-dike scenario (Figure 7l).Additionally, more infiltration occurred for the dike scenario compared to the no-dike scenario after 4 years (Figures 7g,7h,7k,and 7l).Lastly, the salinity following the 2-m surge with the dike removed (Figure 7k) is higher compared to the same surge height with the dike present (Figure 7g), as more salt water inundates the land when the dike is removed.However, the infiltration depth after 1-year post-surge is nearly the same for both the dike and no-dike scenarios (Figures 7g and 7k).
To investigate the influence of tidal pumping on aquifer flushing rates, the same 2-m surge simulations were run post-surge with a mean sea-level boundary (i.e., runs 12 and 13, Table 2) to compare to the scenarios with tides imposed (runs 10 and 11, Table 2) following the same surge (Figures 7i, 7j, 7m, and 7n).The results for the 2-m surge with (Figures 7i and 7j) and without (Figures 7m and 7n) the dike (and without the tide) both indicate the inundated seawater migrates uniformly downward, except near the dike However, for the dike scenario without the tides (Figure 7j), there is deeper infiltration with more diluted seaward that is less than the 1% threshold.In contrast, for the scenario with the dike removed and without the tides (Figure 7m), little dilution and infiltration has occurred.1).(a) Pressure head distribution, (b) salt concentration (normalized to seawater) distribution and wedge location limit based on monitoring well electrical conductivity, and (c) the simulated and observed monitoring well groundwater levels (Figure 2c).Model domain is vertically exaggerated to visualize changes.
Aquifer SWI and subsequent recovery results with (Figures 7g and 7h) and without (Figures 7i and 7j) the tidal pumping included are similar; however, density-driven salt fingers (e.g., Post & Houben, 2017) began to develop in the tidal case (Figures 7g and 7h), whereas infiltration was more uniform for the no-tide case (Figures 7i and 7j).These patterns would likely persist if simulations were run out for longer than 4 years, indicating that the mega tides may accelerate aquifer recovery, in agreement with other studies that have shown even small-amplitude tides increase solute mixing and spreading (e.g., Pool et al., 2014).The increase in saltwater volume was 10.4% and −39.0% for the tide and no tide scenarios, respectively (Figures 7g and 7i) after 1 year post-surge.The "no tide" scenarios have more freshwater than the initial conditions because the USP disappears without tidal pumping.
For the results with the dike removed and with (Figure 7k) and without (Figure 7m) tides, the simulated infiltration pattern (both depth and concentration) is generally similar among scenarios.However, the overall saltwater volume increase was 22.1% and −31.2% for the no dike and no tide/dike scenarios, respectively (Figures 7k  and 7m) after 1 year post-surge, and 44.5% and −57.3% for the no dike and no tide/dike scenarios, respectively (Figures 7l and 7n) after 4 years post-surge.1), and (c) 1-year with the tides back on after SLR (run 5, Table 1).Solid white lines indicate the 1% isocontour location at the initial condition of the SLR runs, and black solid lines indicate post-SLR 1% isocontour location.Results outputted at mid-tide stage.

Discussion of Model Results
With a high connectivity between the basin and aquifer, any ocean perturbations such as SLR (Werner et al., 2013), intensifying storm surges (Bevacqua et al., 2020;Bhatia et al., 2019;Tebaldi et al., 2012) or storm tides (Desplanque & Mossman, 1999;McLaughlin et al., 2022;Sebastian et al., 2014), or anthropogenic impacts from increased pumping rates (Almar et al., 2021;Ferguson & Gleeson, 2012;Klassen & Allen, 2017) could all profoundly impact the availability and use of freshwater resources in these coastal areas.Results from the SLR simulations reveal that the USP stability is highly sensitive to maximum sea level during flood tides, even for very small increases in the maximum high tide.These findings relate to the USP sensitivity to seasonal external forcing (i.e., changes in freshwater discharge), which can dynamically transform the USP and shift the conditions from stable to unstable (Fang, Zheng, Guo, et al., 2022).Similar to how salt wedge migration occurs due to SLR (e.g., Werner et al., 2013), these results highlight that the USP also migrates landward in response to an elevated sea level (Figure 6c).In the context of wellfield salinization, upconing or saltwedge migration have received by far the most research attention.Our results suggest that USP migration may also threaten drinking water resources and force well relocation, at least in areas with high tidal ranges and thus large USPs.
In general, the storm surge simulations were conservative given the simplified repeating single-harmonic tidal boundary condition employed (Section 4.3.1).A much lower surge height (e.g., <1 m) would be required for dike overtopping, if the surge was coincident with spring-high tides (bimonthly cycle, McLaughlin et al., 2022) and even lower (e.g., <0.5 m) if coincident with the saros tide (18.03 years cycle, Desplanque & Mossman, 1999;Gonneea et al., 2013).Additionally, the surge heights used in this study were relatively small compared to other storm surge and overtopping case studies, as other studies were often based on tsunamis or tropical cyclones (Cardenas et al., 2015;Illangasekare et al., 2006), or used several recurring surges of increasing magnitude and frequency (Goldenberg et al., 2001;Paldor & Michael, 2021;Tebaldi et al., 2012).In contrast, this study only simulated a single, short-term surge event for different scenarios (with/without tides and dikes), as the Bay of Fundy region does not typically experience large surge or high-wave events, and this region is more vulnerable to storm tides when high tides and surges are coincident (McLaughlin et al., 2022).Removing the dike only exacerbates the risk of SWI from surges extending further inland.Once the saline water recharges an aquifer it may be pumped landward if it enters a wellfield capture zone.These potential subsurface impacts of intentional dike breaching have generally been overlooked in coastal management given the focus on using nature-based solutions to promote other (surface) ecosystem services.
Shallow soil salinization from dike overtopping or removal could result in the aquifer being contaminated with saline water for decades, even for only short surge durations (Chang et al., 2011;Morgan et al., 2015;Watson et al., 2010;Yang et al., 2018).Many factors influence aquifer recovery (flushing) following vertical SWI, including density driven phenomena, groundwater recharge, and land-sea hydraulic gradients (Cantelon et al., 2022).
For example, salt fingering may create localized circulation cells that influence groundwater flow and salinization dynamics (Xiao, Li, et al., 2019;Xiao, Wang, et al., 2019).Thus, while most studies of vertical SWI due to coastal flooding have been conducted for rapidly flushed beach aquifers, heterogeneous groundwater systems with surficial lower-permeability units, as in the present study, may take much longer to salinize and to recover (e.g., Yang et al., 2013).Additionally, erosion of coastal barriers from storm surges, while not considered in this study, can also impact SWI and aquifer recovery dynamics (Cantelon et al., 2022;Elsayed & Oumeraci, 2018c).Storm surges also threaten certain agricultural species and fields located within kilometers of the coast (Nicholls et al., 2021;Pitman & Läuchli, 2002;Schieder et al., 2018).Due to crop intolerance to salinity and/or increased duration of soil saturation during the growing season, producers along the coast will likely change the kinds of crops they grow, or abandon the farmland for coastal marshlands (Butcher et al., 2018;Fagherazzi et al., 2019;Gedan & Fernández-Pascual, 2019;Guimond & Michael, 2020).
While there are economic and cultural implications associated with abandoning agricultural activities immediately behind the dikes, there are also some potential benefits which include rewilding coastal ecosystems (i.e., wildlife habitat) and carbon sequestration (Hopkinson et al., 2012).Tidal wetland expansion through managed realignment of dikes can provide flood protection and wave attenuation (MacDonald et al., 2017), but may also deleteriously impact crop production in former dikelands and may salinize former groundwater drinking supplies in the long term (Tackley et al., 2023).Dike management decisions are multi-faceted and involve many stakeholders (Sherren et al., 2021), and no single solution can protect or maximize all ecosystem services at play.Other factors not considered in this study, including changing groundwater recharge, pumping rates, or land subsidence, may further complicate future groundwater dynamics in these dikelands (e.g., Ferguson & Gleeson, 2012;Werner et al., 2013).However, findings from this study show that dike management techniques (e.g., dike removal) directly impact coastal freshwater and agricultural soils.Therefore, decisionmakers should explicitly consider how dike management techniques will influence soils and coastal aquifers used for drinking water and irrigation and communicate these implications to dikeland stakeholders.

Model Simplifications and Recommendations
Given the complexity of the modeling activity undertaken (mega-tidal forcing with variably-saturated, variable-density flow), simulations were run using a super-computer network cluster (Siku ACENET) but still required run times exceeding several weeks.This ultimately limited opportunities for adding additional complexities or sensitivity analyses, although the calibration process described in Section 4.3.1 did entail a detailed investigation of the impacts of hydraulic and storage properties on the salinity and head propagation replicated in the model.When paired with the new field data, this modeling study provides a critical baseline for future studies, and the following recommendations are provided to expand upon this study, particularly as computational power continues to increase.
To improve understanding of the site conditions at the study site or in other coastal aquifers, more extensive geological (e.g., deep boreholes with logs), geophysical (e.g., downwell or offshore to better map the USP or the offshore freshwater-saltwater interface), and hydrogeological data (e.g., spatial and temporal hydraulic head and EC data) could be collected to help constrain model parameterization.While no active pumping was conducted in this study, future studies could consider how USP migration in settings with high tidal ranges may be influenced by groundwater extraction.To improve model calibration, the unsaturated zone in the dikeland till and the salt marsh in front of the dike could be parameterized using field data.Future studies could also be run out longer post-surge, to assess long-term flushing rates, and to consider the impacts of recurring flooding (e.g., Paldor & Michael, 2021) or multiple SLR rates based on different climate models or RCPs.Future work might consider different dike locations (managed realignment, Tackley et al., 2023), as the present study only considered the presence or absence of dikes as "bookend" dike conditions).

Summary and Conclusions
This study investigated the potential for SWI in the dikelands along the mega-tidal Bay of Fundy in Nova Scotia, Canada, in response to SLR, storm surge overtopping, and dike removal.Results from the field investigations highlight pronounced tide-driven bay-aquifer interactions as revealed by high-amplitude tidal oscillations in the deep groundwater head over 100 m inland.Additionally, the dynamic intertidal region was captured by observed intra-tidal variability in subsurface electrical resistivity reaching to several 10's of meters.This study, which is the first to consider groundwater-surface water interactions along a mega-tidal coastline, revealed that these overlooked subsurface exchanges are pronounced and may influence shoreline or benthic biogeochemistry.
Results from the numerical model simulations show that the coastal aquifer is vulnerable to increased salinization (e.g., USP expansion and instability as well as vertical SWI) from both rising seas and storm surge events.Upper saline plume dynamics have received less attention than saltwater wedge migration dynamics in the SWI literature, and our results show that migration of these zones may be important to consider as a SWI mechanism along coastlines with high tidal ranges.Vertical SWI can occur during and following storm surges, and the model results indicate that salinization and freshening patterns are influenced by tides as well as dikes.This climate change assessment ultimately highlights new risks that should be considered in dike management and climate change adaptation plans as most dikeland decision making has focused on surface water impacts.Further research and knowledge transfer related to the unseen subsurface influence of tides, dikes, SLR, and storm surges is needed to increase stakeholders' and decisionmakers' awareness and to develop climate-resilient coastlines worldwide.We acknowledge the support of the Natural Sciences and Engineering Research Council (NSERC) ResNet project, funding reference number NSERC NETGP 523374-18.N. LeRoux was supported by a Nova Scotia Graduate Scholarship and by funding from Agriculture and Agri-Food Canada project J-002305-Environmental Change Onehealth Observatory (ECO 2 ) to D. Lapen.HydroGeoSphere modeling technical support was provided by Aquanty.Funding for this project was also provided from the Canada Research Chairs Program and an NSERC Discovery Grant held by B. Kurylyk.Advanced computing resources were provided by the Digital Research Alliance of Canada, the organization responsible for digital research infrastructure in Canada, and ACENET, the regional partner in Atlantic Canada.Lastly, we thank the Wolfville municipality (Director of Engineering & public works and plant operator) for their support in data collection for this study.We appreciate insightful comments from the Associate Editor Larry Murdoch, Marc Walther, and two anonymous reviewers who helped improve the manuscript.

Figure 1 .
Figure 1.(a) Conceptual model of coastal aquifer conditions in a dikeland with a mega-tidal range.The upper saline plume (USP) is delineated.Other panels indicate (b) the influence of mega tides on aquifer salinity and lateral saltwater intrusion (SWI) via (c) sea-level rise and pumping, and (d) vertical SWI from overtopping.

Figure 2 .
Figure 2. Study site map indicating (a) the relative location in Canada, (b) the location within Nova Scotia, and (c) the study site instrumentation near the town of Wolfville.

Figure 3 .
Figure 3. (a) Model domain with surface and subsurface boundary conditions and (b) flow chart of modeling process.Superscript numbers in the flow chart relate to run numbers in Table1.

Figure 4 .
Figure 4. Field data collected in or near the Wolfville dikelands (Figure 2c).(a), (b) Tidal surface water levels (light blue), surface water electrical conductivity (EC) (light orange), groundwater well level (dark blue), and groundwater EC (dark orange), (c) results from the frequency domain analysis of the surface water and groundwater level timeseries, highlighting the dominant higher-frequency tidal constituents (dashed lines), (d) time-lapse bulk resistivity survey collected (e) every hour from ∼low tide to ∼spring-high tide on 25 June 2021 (colors in a, b, and c line up, and colors in d and e are in alignment).

Figure 5 .
Figure 5. HydroGeoSphere model calibration results (run 3, Table1).(a) Pressure head distribution, (b) salt concentration (normalized to seawater) distribution and wedge location limit based on monitoring well electrical conductivity, and (c) the simulated and observed monitoring well groundwater levels (Figure2c).Model domain is vertically exaggerated to visualize changes.

Figure 6 .
Figure 6.Salt distributions as a result of lateral saltwater intrusion from sea-level rise (SLR) after (a) 1 year of SLR with the tides off, (b) 28-years of SLR (i.e., by 2050, run 4, table1), and (c) 1-year with the tides back on after SLR (run 5, Table1).Solid white lines indicate the 1% isocontour location at the initial condition of the SLR runs, and black solid lines indicate post-SLR 1% isocontour location.Results outputted at mid-tide stage.

Figure 7 .
Figure 7. Salt distributions following storm surge events (a) 24 hr after the 1 m surge (run 6), (b) 1 year after the 1 m surge (run 9), (c) 24 hr after the 2 m surge (run 7), (d) 1 year after the 2 m surge (run 10), (e) 24 hr after the 2 m surge with the dike removed (run 8), (f) 1 year after the 2 m surge with the dike removed (run 11), (g) zoomed in panel from (d), (h) zoomed in panel from (d) ran for 4 years post-surge, (i) zoomed in panel from (d) but run without tides, (j) zoomed in panel from (d) but run without tides for 4 years following the 2 m surge, (k) zoomed in panel from (f), (l) zoomed in panel from (f) ran for 4 years post-surge, (m) zoomed in panel from (f) but run without tides, and (n) zoomed in panel from (f) but run without tides for 4 years following the 2 m surge with the dike removed.Solid white lines indicate calibration 1% isocontour location (initial outside edge of upper saline plume), and black solid lines indicate post-surge 1% isocontour location.Results outputted at mid-tide stage.

Table 2
Properties Used in HydroGeoSphere Model Simulations