Field-scale modelling reveals dynamic groundwater flow and transport patterns in a high-energy subterranean estuary

Subterranean estuaries (STEs) are biogeochemical reactors modifying the chemistry of salt- and freshwater as they flow through the subsurface sediments. Boundary conditions such as tides, waves, beach morphology, seasonal meteoric groundwater recharge and storm events control endmember mixing and residence time distributions within STEs. These in turn affect biogeochemical reactions and thus elemental fluxes discharging to the ocean via submarine groundwater discharge (SGD). Especially at high-energy beaches exposed to high tidal ranges and high wave energy, boundary conditions are very dynamic and likely imprint on groundwater flow and reactive transport within the STEs. A quantitative understanding of mixing processes and residence time distributions is necessary in order to adequately describe biogeochemical processes and can be achieved with the help of numerical modelling. Yet, transient field-scale modelling approaches calibrated to comprehensive observational data sets are still lacking, in particular for real-world high-energy STEs. In the present study, for the first time a density-dependent groundwater flow and transport model was developed and calibrated for a high-energy beach. The north beach of the barrier island Spiekeroog, northern Germany, thereby served as an example field site exposed to high-energy characteristic boundary conditions. The model was calibrated to a 1.5-year extensive dataset of groundwater heads, salinities, temperatures and ³H $/$ He groundwater ages at various shore-perpendicular locations along the beach at depths down to 24 m below ground surface. The calibrated model is able to replicate the principal behaviour of the highly transient system and enabled the identification of hot spots of high temporal variability in the investigated state-variables. The dynamics in salinity are most intense at the in- and exfiltration locations of the tide-induced recirculating seawater. The groundwater age variability was largest seawards of the low tide mark as well as below the deep recirculating seawater cell at around 20–30 m depth near the dunes, where very old freshwater from the islands' freshwater lens mixes with young brackish water from the upper beach. Temperature variations were seasonal and confined to the upper 5–10 m below the beach. Computed saline SGD water fluxes varied considerable on daily and spring-neap time scales, as well as on the longer term, i.e., monthly to yearly time scales. The rather gradual, longer-term changes in flux appear to be mainly controlled by changes in spatial variability of the beach slope. The simulated groundwater age of the fresh SGD component varied between 4 and 25 years, and predominantly depended on the magnitude of saline SGD flux. Overall, the model provided important insights into the dynamics of the flow and transport processes.

1 Introduction

Subterranean estuaries (STEs), or coastal beach aquifers, control the fluxes of carbon, nitrogen, phosphorous and other elements delivered to the global ocean via submarine groundwater discharge (SGD), and play an important role in global biogeochemical cycles (Santos et al., 2021; Wilson et al., 2024). Quantification of these elemental fluxes, however, is difficult, because both (i) fresh and saline SGD water fluxes are not easy to measure under field scale conditions (e.g., Burnett et al., 2006; Grünenbaum et al., 2020a) and (ii) complex hydrogeological and hydrobiogeochemical processes, and thus reactive transport within the STE, modify solute concentrations prior to discharge (Anwar et al., 2014; McAllister et al., 2015; Kim et al., 2017; Heiss et al., 2017; Kim et al., 2019, Cogswell and Heiss, 2021). These processes are affected by the physical (permeability and porosity) and biogeochemical (e.g., abundance of reactive minerals, mineral surfaces and particulate organic matter) properties of the beach aquifer (Kim et al., 2017; Heiss et al., 2017, 2020; Seibert et al., 2024). Another impact on the reactive transport system are highly dynamic hydrological boundary conditions such as tides, waves, storm surges and seasonal terrestrial groundwater recharge (Santos et al., 2012; Robinson et al., 2018; Greskowiak at al., 2023). Transient hydrogeochemical boundary conditions, e.g., water temperatures, concentrations and lability of dissolved organic matter, and other reactants which are present in the fresh- and saltwater end members (Ahrens et al., 2020; Greskowiak et al., 2023; Seibert et al., 2025) affect these processes even further. Subsurface water residence times, flow paths, and dispersive mixing variable in space and time control the extent to which hydrobiogeochemical transformations can occur (Robinson et al., 2009; Anwar et al., 2014; Heiss et al., 2017). Thus, thorough understanding of groundwater flow and solute transport is critical for the identification and quantification of reactive processes within the STE. Numerical modelling thereby aids to improve the understanding, as non-intuitive system behaviour stemming from non-linearities (e.g., density effects) and/or superimposition of dynamic boundary effects can be captured, disentangled and analysed (Meyer et al., 2025a). Therefore, in STE research, numerical groundwater flow and transport modelling has already been applied in many studies.

Lebbe (1981) was likely the first to observe and model the saltwater-freshwater distribution within a tide-affected beach aquifer at the De Panne field site (Belgium). Observations were made by up to 30 m deep electrical resistivity profiles measured at 5 shore-perpendicular transects, each consisting of 4 to 6 boreholes from dunes to the low tide mark, as well as piezometric head data at these locations. They revealed that in the intertidal zone, the beach aquifer was characterized by a seawater lens overtopping freshwater with a density-induced saltwater wedge below. Decades later, the overtopping seawater lens and underlying freshwater were termed upper saline recirculation cell or upper saline plume (USP, Robinson et al., 2006) and freshwater discharge tube (FDT, Robinson et al., 2006), respectively. Lebbe aimed to describe and interpret these observations with the help of numerical modelling. His model was a finite-difference, density-dependent, steady-state, 2D vertical cross-sectional model. It employed a streamline approach that iteratively delineated the zones of freshwater and seawater, and a tide-averaged ocean boundary condition depending on intertidal topography. The model could replicate the principal saltwater-freshwater distribution at all transects studied, and deduce the flow patterns in the different compartments of the STE. Later, other field-site investigations and accompanied density-dependent flow and transport modelling (e.g., Lebbe, 1999; Robinson and Gallagher, 1999; Robinson et al., 2006; Vandenbohede and Lebbe, 2007) confirmed the principal configuration of USP, FDT and salt water wedge of tide affected STEs. However, numerical simulations of Evans and Wilson (2016) showed that USP formation strongly depends on the prevailing hydraulic conductivity and beach slope, and may be less common than previously thought. Other modelling studies revealed the importance of subsurface heterogeneity (Geng et al., 2020a, b; Olorunsaye and Heiss, 2024), spring-neap tide cycles (e.g. Robinson et al., 2007; Abarca et al., 2013; Heiss and Michael, 2014), intensified wave conditions by offshore storm events (e.g., Robinson et al., 2014) and inland seasonal meteoric groundwater recharge (Michael et al., 2005) on the formation of the USP and its dynamics, as well as on groundwater flow patterns, SGD fluxes and subsurface travel time distributions. Also, storm events (e.g., Holt et al., 2019) and time-varying intertidal morphology (e.g., Heiss and Michael, 2014; Greskowiak and Massmann, 2021) can considerably affect the subsurface flow and transport dynamics in the STE (Meyer et al., 2025a). Comprehensive reviews of the different impacts on subsurface flow and transport dynamics in STEs are given by Robinson et al. (2018) and Geng et al. (2021).

The vast majority of the existing model applications in the current STE literature are of generic or semi-generic nature, considering idealized aquifer-settings and/or boundary conditions in order to investigate a specific aspect or an individual (sub)process of the system. Nevertheless, studies exist that applied a realistic field-scale model and calibrated it, or at least compared it to observations from a real field site (Robinson et al., 2006; Vandenbohede and Lebbe, 2007; Abarca et al., 2013; Heiss and Michael, 2014; Evans and Wilson, 2017; Zhang et al., 2017; Beck et al., 2017; Grünenbaum et al., 2020b; LeRoux et al., 2023). But of those, many investigations were carried out at protected beaches with low or moderate energy conditions, while modelling studies of high-energy beach aquifers are scarce. The existing studies considered only idealized boundary conditions and dynamic equilibrium of groundwater flow and salinity distribution (Vandenbohede and Lebbe, 2007; Beck et al., 2017; Grünenbaum et al., 2020b; LeRoux et al., 2023). Typically, these studies suffer from limited data availability for site characterization and model verification. This is because the intertidal zone, especially at a high-energy site, is hard to instrument due to the harsh environmental conditions (Anschutz et al., 2009; Charbonnier et al., 2013; Massmann et al., 2023). The lack of comprehensive modelling applications for high-energy sites may have also been caused by the difficulty to numerically deal with highly transient and irregular forcing at the model boundaries. In combination with memory effects lasting several months (e.g., Robinson et al., 2014) and the non-linearity of the system, this requires very large simulation and computation times, making model calibration cumbersome (e.g., Holt et al., 2019).

The aims of the present study are to (i) provide a numerical groundwater flow and transport modelling analysis for a high-energy beach site that can explain real-world field observations and, thereby, (ii) disentangle the complex system behaviour and reconstruct the physical processes in the beach aquifer as a consequence of various superimposing hydro(geo)logical forces in a real high-energy beach setting, (iii) identify potential hotspots of intense gradients and variability, and (iv) estimate the fresh and saline SGD fluxes and respective mean subsurface residence times, as well as to illustrate their transience and the reasons for it. The high-energy beach of the barrier island Spiekeroog in Northern Germany (Fig. 1), which is presently intensely studied (Massmann et al., 2023), will serve as an example field site for that purpose. An extensive and unique tracer dataset, consisting of continuous groundwater head, salinity and temperature measurements, as well as occasional groundwater age dating (Meyer et al., 2025b; Skibbe et al., 2024; Reckhardt et al., 2024) is used for model calibration and validation. The model presents the vital basis for subsequent reactive transport modelling and quantification of biogeochemical transformation rates at the Spiekeroog case study site.

Figure 1 Location of the instrumented field site on Spiekeroog Island. The red line indicates the model transect. Dark grey area on the island highlight dune areas. The red circles indicate the positions of the multi-level groundwater monitoring wells ML1, ML2 and ML3. The green circle represents the geophysical saltwater monitoring system SAMOS. The yellow circles indicate the positions of the Direct-Push (DP) sampling points. The aerial picture is from a drone-survey on 26 June 2022. The maps were adapted from Massmann et al. (2023) under the open access Creative Commons Attribution License (CC BY).

2 Methods

For modelling the density-dependent groundwater flow and transport processes the USGS code SEAWAT (Langevin et al., 2008) was used. SEAWAT couples the groundwater flow code MODFLOW (Harbaugh, 2005) and transport code MT3DMS (Zheng and Wang, 1999) to account for density effects on groundwater flow, and it has been proven as suitable tool for the simulation of flow and transport in STEs (e.g., Robinson et al., 2006; Heiss and Michael, 2014). The governing equations are given in the Supporting Information. Advection-dispersion was solved with the implemented Modified Methods of Characteristics (MMOC) scheme, which produces very low numerical artefacts independent of the Peclet number and thus spatial discretization and dispersivity (Zheng and Bennett, 2002).

Model parameterization was based on measured aquifer parameters and long-term hydrological data of groundwater recharge and oceanographic forcing. The simulation time was 21 years with a 19.5-year spin-up period starting on the 1 January 2003. The last 1.5-years served as calibration period, ending in January 2024. The simulation included density-dependent groundwater flow, advection and dispersion of salinity, temperature via heat transport, as well as groundwater age. The latter is a measure of the time that a water parcel resides in groundwater after entering the aquifer, either via meteoric groundwater recharge or recirculating seawater.

For model calibration various datasets were available. At three multi-level groundwater monitoring wells that are aligned along a transect from dunes to low tide water mark (Fig. 1) and screened at 6, 12, 18, 24 m depth, time-series of hydraulic head, salinities and temperatures were acquired (Meyer et al., 2025b, c). Further, a geophysical electrode chain (saltwater monitoring system SAMOS, Ronczka et al., 2020) near the high tide mark recorded subsurface temperatures every 30 min at a spatial resolution of 2 m down to 20 m depth (Skibbe et al., 2024). Moreover, groundwater from all observation wells was dated twice within this period using ³H $/$ He age dating (Meyer et al., 2025b). In addition, high-resolution data of salinity at locations between the multi-level observation wells was obtained from Direct Push groundwater sampling every 12 weeks (Reckhardt et al., 2024).

2.1 Model domain and spatial discretization

The model domain represents a 2D vertical cross-section for the sandy beach aquifer at the northern beach on Spiekeroog Island (Figs. 1 and S1 in the Supplement), oriented parallel to the approximate flow direction from the interior of the island to the sea. In the horizontal direction the model is 1450 m long, extending from the estimated island's groundwater divide in the south (Hähnel et al., 2023) to approximately 200 m seaward of the topography-dependent position of the mean low tide mark in the North (Figs. 2 and S1). The horizontal discretization varies between 150 and 20 m in the first 300 m from the groundwater divide, between 10 and 2 m towards the dune base at 850 m from the divide, and is set to 2 m between dune base and the northern model boundary. The maximum depth of the model domain is defined at 40 m below the mean sea level, where a continuous clay forms the basis of the aquifer (Röper et al., 2012). In vertical direction, the model consists of 80 layers. Except for the uppermost layer, the vertical discretization ranges from 0.51 m at the southern to 0.42 m at the northern boundary. From the groundwater divide towards the dune base, the top of the uppermost layer follows the topography. From the dune base onwards, it decreases with a constant slope from 3.5 m above mean sea level (m a.s.l.) to −6 m a.s.l. at the northern boundary. The vertical extent of the uppermost layer had to be defined as such that it cannot fall dry during the simulation in case the groundwater-level drops in response to transient changes at the sea-level boundary. This was done to avoid the SEAWAT-specific numerical difficulties that sometimes occur during frequent drying and re-wetting of model layers in transient phreatic aquifer simulations. The thickness of the uppermost layer is 1 m, where the top of the layer is equal or less than 1 m a.s.l. Where the top of this layer was above 1 m a.s.l., the bottom stays constant at 0 m a.s.l., and as a consequence, the layer thickness continuously increases towards the dune base, where it reaches 3.5 m. The saturated thickness of the first layer, however, typically ranges only between 1 and 1.5 m (depending on season and tides), but can occasionally increase to the maximum layer thickness in the upper beach in response to storm events and associated strong rise in groundwater levels (see observed hydraulic head at multi-level observation well ML1 at 6 m depth, Fig. S2 in the Supplement). Groundwater flow is simulated as unconfined where the hydraulic head is below the top of the first layer. If it is above, confined conditions are applied.

Figure 2 Model-domain including the calibrated horizontal hydraulic conductivity field as well as the locations of pilot points for model calibration (black squares), multi-level observation wells ML1, ML2 and ML3 (white dots) and the geophysical electrode chain SAMOS (red dots). Horizontal and vertical grey lines show the model grid. The geological structure was adopted from Massmann et al. (2023). Direct Push sampling locations are not shown here.

2.2 Boundary and initial conditions

The bottom boundary and the two vertical boundaries at both sides of the model domain are defined as no-flow boundaries (Fig. S1). The top boundary received either monthly varying meteoric groundwater recharge (Fig. S1), with an annual average of 425 mm a⁻¹ (Hajati et al., 2022), via a 2nd type boundary condition (details on the recharge calculation are given in the Supporting Information, including the applied recharge time series in Fig. S3), or a defined hydraulic head via a 3rd type boundary condition for nodes occupied with seawater.

At the sea boundary, a tide-averaged hydraulic head was defined, as tide-resolved computations for a simulation time of 20 years are computationally too demanding when aiming for model calibration, even with a 2D cross-sectional model. The applied hydraulic heads were a function of local topographic height, tidal signal and mean significant wave height, based on the tide-averaging approach of Nuttle (1991). The approach of Nuttle has already successfully been used in various adaptations to simulate groundwater flow and transport in tide-affected beach aquifers in a computationally efficient way (Vandenbohede and Lebbe, 2007; Pauw et al., 2014; Grünenbaum et al., 2020b; Greskowiak and Massmann, 2021; Hähnel et al., 2023; Seibert et al., 2024). In order to reflect the importance of the wave setup on the STE flow system (e.g. Xin et al., 2010, 2014; Robinson et al., 2014), the approach of Nuttle (1991) was extended by adding the wave setup using an empirical formulation by Nielsen (2009), which has successfully been used before to describe wave-induced groundwater flow below beaches (e.g., Malott et al., 2016; Wu et al., 2017). Details on how this was computed is given in the Supplement. Note that during the receding tide, the groundwater table around the high water mark generally starts to decline due to drainage, which, in a sandy beach aquifer, can be a few decimeters (Vandenbohede and Lebbe, 2007). In Spiekeroog, the average decline over one tidal cycle was found to be about 20 cm (Grünenbaum et al., 2020b). To account for this, the calculated tide-average hydraulic heads at the sea boundary condition were further modified after Grünenbaum et al. (2020b) and Greskowiak and Massmann (2021). Therein it was assumed that at the mean high water mark (which is 1.35 m a.s.l.) and above, the average decline in hydraulic head is 20 cm, while in seaward direction it linearly decreases to zero when reaching the mean sea level mark.

Data on tides and offshore significant wave height was provided by the German governmental service of maritime navigation and the seas (BSH). For the first 10 years of simulation time, i.e., from 1 January 2003 until 31 December 2012, an average semi-diurnal tide with an amplitude of A=1.35 m and an average offshore significant wave height of H_sig=0.58 m, as well as a constant beach morphology from May 2014 (Beck et al., 2017) was used to calculate the average hydraulic heads at the aquifer-sea boundary. From 2013 onwards, tide data at one-minute resolution was used for the calculation of phase-averaged hydraulic heads over a tidal period, i.e., from low water time to the low water time of the next tide, according to the scheme described above. High tide and low tide data is presented in Fig. S4 in the Suppleemnt. Observed 1 min data on offshore significant wave height was available from 12 August 2020 onwards. Before that date, 3-hourly WaveWatch III hindcast model data (The WAVEWATCH III^® Development Group, 2019) was used. The employed wave data set is presented in Fig. S5 in the Supplement. The beach morphology from the dune base to the low tide mark was irregularly measured during low tide with a differential GPS (vertical accuracy: ± 1.5 cm) between May 2014 and March 2022 (Fig. S5 in the Supplement). Since April 2022, a beach profile was measured every 6 weeks. Between the low tide mark and the northern model boundary, the topography was linearly interpolated in space. From 2013 onwards, the tide-averaged hydraulic heads at the aquifer-sea boundary were computed for each tidal cycle in so-called stress periods, based on the time series of tide, significant wave height and beach topography. For each stress period, the beach topography was assumed constant and its elevation at each location received a value from linear interpolation in time (Fig. S6). Each stress period was further subdivided into three time steps to accommodate for the transient response of groundwater flow due to the changing hydraulic heads at the aquifer-sea boundary. Note that a very strong storm flood eroded large parts of the beach between 15 and 17 February 2022, thereby shrinking the width of the intertidal zone by around 100 m (Fig. S6). Unfortunately, topography measurements before and after this event were only available from July 2021 and April 2022, respectively, as this was the time before the intense sampling. In the model, it was assumed that the topographic changes in-between these dates occurred instantaneously during the time of the storm event.

Transport of salinity, groundwater age and temperature were simulated with the model. Therein, the groundwater recharge and aquifer-sea boundaries were assigned a non-dispersive flux boundary condition, and thus mass is transported in or out of the aquifer only via advective flow at these boundaries. At the groundwater recharge boundary, a salinity of zero, a groundwater age of zero and a water temperature following a sine function with a minimum of 4 °C at the end of January and a maximum of 21 °C at the end of July was applied. The same was done for the aquifer-sea boundary at locations were water entered the model domain, except that a salinity of 32 (PSU, Practical Salinity Unit) was applied. Both recharge and aquifer-sea boundaries were assigned the same seasonal temperature variation based on observed shallow seawater temperatures (Reuter et al., 2009), since temperatures of the groundwater recharge were unknown. Due to heat diffusion, the modelled temperature in the freshwater body was quickly damped to around 10–11 °C by 2–3 m below the groundwater level which is consistent with the observed constant groundwater temperature of the island (Röper et al., 2014; Seibert et al, 2018).

The initial conditions for hydraulic head and temperature were set to 1 m and 10 °C in the entire model domain, respectively. The initial salinity was set to zero in the freshwater lens and the beach area where the topography was higher than the mean spring high tide mark of 1.78 m a.s.l. In the remaining model domain, initial salinity was set to 32. The assigned initial groundwater age linearly increased with depth from 4 to 51 years according to ³H $/$ He age dating within the island's freshwater lens (Röper et al., 2012; Seibert et al., 2018).

2.3 Parameterization and model calibration

Model parameters of the groundwater flow model were the horizontal and vertical hydraulic conductivities K_H and K_V, the specific yield S_y, specific storage S_s, as well as the so-called hydraulic conductance C_h (m d⁻¹), which describes the degree of hydraulic connection between aquifer and sea via a 3rd type boundary condition. The latter was chosen to simulate the response of the aquifer to the external forcing of the sea, because it has been proven as numerically very stable in tidal beach groundwater simulations (Mulligan et al., 2011). The hydraulic conductance was defined as $C_{\mathrm{h}}=c_{\mathrm{scale}}\cdot (\mathrm{\Delta }x\cdot K_{\mathrm{V}})\cdot {\mathit{\delta }}^{-\mathrm{1}}$ (Mulligan et al., 2011; Hähnel et al., 2023), where δ is the vertical distance between ground surface and computational node, and Δx the horizontal width of the boundary cell. However, the computation of C_h is not straightforward, as the flow direction is not known a priori (Mulligan et al., 2011). Thus, C_h was subject to calibration via c_scale, which was allowed take a value between 0.3 and 3.0. For transport simulations, the model parameters were the effective porosity n_e, as well as longitudinal and vertical transverse dispersivities α_L and α_VT, respectively. All of the above model parameters were subject to calibration, except the specific storage, which was fixed to $S_{\mathrm{s}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m⁻¹. In addition, for the simulation of temperature, a fixed thermal retardation factor of R=2.0 and a bulk thermal diffusivity of D=0.2 m² d⁻¹ were assigned, in agreement with other groundwater flow and transport studies in sandy aquifers (e.g., Prommer and Stuyfzand, 2005; Greskowiak et al., 2006; Thorne et al., 2006; Sharma et al., 2012; Henzler et al., 2016). The groundwater age simulation was performed after the direct groundwater age approach of Goode (1996), in which a solute component is assigned a zero-order production rate k=1 T⁻¹, where T is the time unit. The density of water ρ was only made dependent on salinity (see governing equations in the Supporting Information), which is a generally accepted simplification of the equation of state for density-dependent groundwater flow in seawater affected freshwater aquifers (e.g., Thorne et al., 2006). Note that preliminary simulations showed that the temperature-dependence of both water density and viscosity on the overall flow patterns in this system was insignificant due to the rapid dampening of the temperature signal.

The subsurface at the study site consists of three more or less horizontally layered geological units, i.e., beach sands from the ground surface to approximately −3 m a.s.l., tidal flat deposits until around −15 m a.s.l., and glacio-fluvial sediments at greater depth (Massmann et al., 2023; compare Fig. 2). Slug and Bail tests in the observation wells, infiltration tests at the beach surface, as well as grain size analyses and Darcy experiments with sediments from cores retrieved during the installation of the monitoring wells gave first estimates on horizontal and vertical hydraulic conductivities (Meyer et al., 2025b). The horizontal hydraulic conductivities of beach sand K_{H_B}, tidal flat K_{H_TF}, glacio-fluvial deposits K_{H_GF} were in the range of 15–45 m d⁻¹, 4–40 and 30–80 m d⁻¹, respectively. The vertical anisotropies An $=K_{\mathrm{H}}/K_{\mathrm{V}}$ varied largely between 1 and 4, and total porosities were found in the range of 25 %–40 %. All model parameters subject to calibration are given in Table A1 in the Appendices, including initial guesses and parameter bounds for calibration, as well as the final calibrated values.

Model calibration was done against observed 1.5-year long time series of salinities, ³H $/$ He groundwater ages, temperatures measured at all monitoring wells as well as temperatures recorded at the vertical geo-electrical chain SAMOS (Figs. 3–5) and hydraulic heads recorded at ML1 at 6 m below ground surface (Fig. S2). Additionally, salinity data from Direct Push sampling (Reckhardt et al., 2024) at the transect was used for model verification.

Automatic parameter estimation was carried out with a particle swarm optimization (PSO) approach (Siade et al., 2019), which is implemented in the model-independent parameter estimation tool PEST$++$ (White at al., 2020). It was the preferred approach because compared to gradient-based optimization schemes, it is less prone to getting stuck in local minima of the objective function, which is a common problem for highly non-linear models and would lead to a premature end of a calibration run (e.g., Siade et al., 2021). A PSO parameter estimation run required approximately 10 000 model simulations and took around 10 d on a High-Performance-Computing Cluster, with running 100 simulations in parallel. A single model simulation required about 2–3 h (depending on parameter values) on an AMD Genoa EPYC 9554 @ 3.1 GHz CPU core. For the optimization, observation weights were applied for 5 individual observation groups, namely (a) hydraulic heads, (b) salinity, (c) temperature at the multi-level monitoring wells and (d) temperature at SAMOS and (e) groundwater ages. In the objective function, the weights are usually set simply as the inverse of the measurement standard deviation. However, in a multi-objective optimization framework, where different types of state-variables with an uneven number of data points in each dataset are used for calibration, weighting is not straight forward anymore (Dai and Samper, 2004): The weights need to balance the uneven amount of data points in the individual datasets of different state-variables and their often vastly different numerical values. In addition, their contribution to the objective function, i.e., the weighted sum of squared residuals, change in the course of the optimization. Thus, an iterative approach is needed (Dai and Samper, 2004). In this work, preliminary optimization runs were carried out to find optimal weights for the different observation groups via trial and error. The final weights are presented in the Supplement. The convergence of the PSO optimization is also presented in the Supplement.

The uncertainty in calculated annual groundwater recharge rate was reflected in the calibration procedure by allowing to take a value of 80 %–110 % of its original value of q_R=425 mm a⁻¹. Elevated beach recharge was assumed as vegetation on the beach is virtually absent and the beach sands are very permeable. Consequently, on the beach, an average annual recharge rate of 600 mm a⁻¹ with a summer recharge of around 230 mm a⁻¹ was defined (see Fig. S2), which was subject to calibration bounded between 80 % and 110 % of its original value. Further, the infiltration rate of seawater into the upper beach sediments during storm events (i.e., when the seawater level is higher than approximately 2.2 m a.s.l.) is very uncertain due to the effects of variably saturated flow and clogging by entrapped air, and an adequate description would require a multi-phase approach of the process. Preliminary calibration runs suggested that the applied 3rd type boundary condition could underestimate the seawater infiltration rate on the upper beach. An adjustment of the hydraulic conductance C_b did not alleviate this behavior. Therefore an additional seawater infiltration q_sea was implemented via a 2nd type boundary condition between actual sea level and mean spring tide high tide mark if the sea level was higher than mean spring high tide. The seawater infiltration q_sea was subject to calibration with a maximum rate of 0.25 m d⁻¹.

Preliminary model calibration runs revealed that a certain degree of subsurface heterogeneity with respect to hydraulic conductivity within the three geological units was needed for a reasonable replication of the observed system behaviour. Heterogeneity was realised via the pilot points approach (Doherty et al., 2010). Therein hydraulic conductivities were defined at a distributed set of defined locations - the pilot points - within each geological unit. These serve as reference points where the hydraulic conductivities were adjusted and used to assign a spatially non-uniform hydraulic conductivity field via Kriging interpolation over the entire model grid at each iteration step during the calibration procedure. Trial and error led to a more or less optimal total number of 37 pilot points for a sufficient degree of heterogeneity (Fig. 2). A lower number of pilot points, and thus less heterogeneity, resulted in a premature end of model calibration with a poorer fit, while a higher number of parameters resulted in an extremely slow convergence rate, not reaching a reasonable fit within a bearable time frame. For the reason of the latter, vertical anisotropies of hydraulic conductivity and effective porosities were assigned as homogenous for each of the geological units, and the longitudinal and vertical transverse dispersivities were set homogenous for the entire model domain.

3 Results and Discussion

3.1 Calibration results

As indicated by the calibration performance measures (Fig. S9 in the Supplement) and visual comparison of the temporal and spatial patterns of the observed and simulated state variables (Figs. 3, 4, and S2), the calibrated model describes the principal observed system behaviour very well. For instance, the effects of storm events in February 2022 and December 2023 on the salinity in ML1–6 m (Fig. 3a) and its decline in between these events were well reproduced. On the other hand, the rather stable salinity at greater depth at this location could not be replicated, as the model overestimated the dynamics (Fig. 3a). The same was true for ML2–24 m, and in addition, simulated salinities were largely too low at this location. However, at ML2–6 m, 12 m and 18 m, the simulated concentrations and their dynamics agree very well with the observations. At ML3–18 m and ML3–24 m, the simulated salinity gradually declined until around July 2023, when it slowly started to rise again, while the observed salinity continuously decreased.

The decreasing trends at ML3-6m and ML3-12m starting in November 2023 and August 2023, respectively, could not be reproduced; here the model underestimated the salinity dynamics. Overall, the simulated salinity agrees well with the data obtained from the Direct Push measurements (Fig. 4a). The vertical gradients were largely reproduced, but horizontal gradients, especially in between ML1 and ML2, were less well replicated. A major discrepancy can be noticed for September 2023, where the observed freshwater-saltwater transition was not matched.

Figure 3 (a) Simulated (red line) and observed data of salinity (blue circles) and (b) temperature (blue line) at multi-level observation wells and (c) temperature (blue line) SAMOS. In panel (a), the black line represents the sea level at high tide and is the sum of high tide and wave setup. In plot B and C, the grey dashed line represents the sea water temperature.

The principal patterns in groundwater age were also captured well by the model (Figs. 4 and 5). However, some deviations between simulated and observed data can be noticed. For instance, near the dune base at ML1–12 m, an elevated groundwater age of around 13 years was found in September 2023, but a groundwater age of only 5 years was computed by the model. As generally groundwater age increases with increasing depth (Vogel, 1967), the observations are rather atypical for reasons unknown. In October 2022, the model age was 11 years at ML1–24 m, while the measured age was only 6 years. Also in ML2–24 m, i.e. in the deeper USP, the model overestimates groundwater age by more than 100 %. Nevertheless, the principally rather young groundwater ages of only a few years within the beach aquifer at this depth as compared to 30-year-old groundwater within the interior of the island (Seibert et al., 2018) could be replicated by the model. Notable deviations also exist at ML3–18 m for September 2023 where the model underestimated the 18-year-old groundwater by around seven years. Groundwater ages in the beach aquifer reflect a mixture of a rather young seawater component and considerably older freshwater components (Fig. 5). It appears that deviations between simulated and observed groundwater age are mainly due to the inability of the model to replicate the observed salinity at these locations, as groundwater ages were overestimated where salinities were underestimated and vice versa (compare Figs. 3a and 5).

Figure 4 Simulated (heat map) and observed (circles) data of (a) salinity, (b) groundwater age and (c) temperature. Dates in the top-right corner show the sampling dates and the corresponding model simulation times. The white arrows represent the groundwater flow directions, and the arrow thickness indicates the pore-water velocity.

The time shifts of the summer and winter peaks in groundwater temperature observed at both the multi-level observation wells and SAMOS were very well reproduced by the model (Fig. 3b and c). However, the dampening by dispersion and heat diffusion was somewhat overestimated with increasing depth, and observed temperature variations at ML1–6 m and ML1–12 m were not well reproduced, the reason being unclear. In the observed data at ML2 it is noticeable that the temperature-based residence times (under consideration of a thermal retardation factor of R=2) are considerably younger than the ³H $/$ He groundwater ages. While between August 2022 and May 2023 the temperature peak shift suggests a seawater residence time of about 1.5 years to ML3–24 m (Fig. 3c), it is almost double according to the ³H $/$ He age dating in November 2022 (Fig. 5). This either suggests that young seawater mixes with older freshwater at these locations, or the transport of the temperature signal by heat-diffusion is more important than advective transport. The latter could be likely due to the fact that most of the infiltrating water is caught by the highly permeable beach sediments leading to a rather shallow high flow zone while less water is transported further downwards (e.g., Fig. 4 and Video S1 in the Supplement). Extracted model-based seawater residence times indeed show that the travel times of the seawater component are larger than those suggested by the temperature peak shifts, and are similar to the observed ³H $/$ He groundwater ages (Fig. 5). This would imply that calculating water residence times from the temperature peak shifts is misleading under these circumstances.

Figure 5 Simulated ages of seawater component (dashed red), freshwater component (dashed blue) and groundwater direct age (solid, varying from blue to red depending on salinity), and observed ³H $/$ He groundwater age (squares). The ages of the sea- and freshwater components were calculated from saltwater-freshwater endmember mixing and the computed direct groundwater age.

The calibration results suggest a heterogeneous subsurface with respect to hydraulic conductivity, which explains some important features in the observed data. For instance, the constantly elevated salinity in ML1–24 m compared to the filter screens of ML1 located above (Figs. 3a and 4a) appears to result from a zone of higher hydraulic conductivity extending from the shallow subsurface between ML1 and ML2 towards the dunes with increasing depth (Fig. 2). Notably, the calibrated hydraulic conductivities, especially for the tidal flat sediments, were either at the lower or upper end of the allowed parameter bounds (Table A1), suggesting that the goodness of fit might have been further improved by allowing wider parameter bounds during calibration. However, given the results of the aquifer tests we refrained to do so.

Further, calibration resulted in up to 20-fold higher vertical anisotropies of the hydraulic conductivity than found in the aquifer tests. Preliminary calibration runs where anisotropy values were forced to not exceed 10 (as suggested from the aquifer tests), resulted in significant salt fingering flow. This generated pronounced two- to three-monthly salinity variations, which were not observed in reality. Similarly, the estimated vertical transverse dispersivity was rather high with α_VT=0.25 m, but employing lower values also led to salt fingering flow. Nevertheless, elevated transverse dispersivities are not unusual in highly transient flow systems. For instance, fluctuating flow directions typically lead to higher apparent transverse dispersivities (Kinzelbach and Ackerer, 1986; Goode and Konikow, 1990). More recently, Rau et al. (2018) reported on elevated apparent dispersion in beach aquifers under wave forcing. Annual meteoric groundwater recharge on the island and beach estimated were 470 mm a⁻¹ and 660 mm a⁻¹, respectively, and hence appear rather high. However, soil-water isotope analysis on the neighbouring barrier island Langeoog revealed similarly high groundwater recharge rates of up to 550 mm a⁻¹ in the dune area (Post et al., 2022).

3.2 Groundwater flow and transport behaviour

The simulation results suggest that groundwater flow and transport within the beach aquifer is dynamic and non-periodic (Fig. 4). This can be seen best in the animation video in the Supplement (Video S1). Supported by the comparable high permeable sediment, flow velocities within the shallow beach sediments are high and change rapidly at a weekly, sometimes daily time scale depending on neap-spring tides and storm events. Salinity varies especially in the area above the mean high tide mark, where mixing with freshwater can occur after spring tides or storm floods. These short-term fluctuations in flow and salinity are superimposed by seasonal effects, stemming mainly from fresh groundwater recharge variations, winter storms and varying beach topography. Seasonal effects can also be traced at greater depth. For instance, the general absence of storm surges during the summer months in combination with enhanced fresh beach recharge leads to dilution of seawater within the USP. The resulting temporary drop in salinity can successively be traced through the USP (ML2 in Fig. 3a, and Video S1), with lag times of a few months (ML2–6 m) to more than one year (ML2–24 m). Superimposed, the sudden shortening of the beach due to the winter storm event in February 2022 led to a sharp but small drop in the otherwise gradual decline in the simulated salinity at ML2 and ML3 during this period (Fig. 3a). The effect of the instantaneous beach erosion event is more visible in the subtidal area seaward of ML3, where remnant seawater is slowly flushed out of the beach aquifer and freshening occurs following the event (Fig. 4a and Video S1). The highest variability of salinity in time also occurred in this area where the fresh discharge tube moves through the subsurface, as well as in the upper beach between ML1 and ML2 (Fig. 6). Here, the salinity gradients are highest and thus flow-induced shifts have the largest impact on salt concentrations at these locations.

Figure 6 Standard deviation of simulated salinity, groundwater age and temperature time series along the modelled cross-section between July 2016 and February 2024.

Groundwater ages below the dunes reflect a typical increase with depth as observed within the interior of the freshwater lens on the island (Seibert et al., 2018) with ages of about 10 years at −10 m a.s.l. and about 30 years at −20 m a.s.l. (Fig. 4b). At the beach, the fresh groundwater flowing seaward is pushed down because of seawater infiltration, dragging older water below the USP. The recirculating seawater within the USP is rather young, around 1–2 years at its deepest location (Figs. 4 and 5). In response to the rapidly changing infiltration and exfiltration patterns and magnitude, the groundwater age distribution in the beach aquifer also varies with time (Figs. 4 and 5, and Video S1). A major effect is visible somewhat upgradient of ML2 from October 2022 onwards, where an area of younger groundwater starts to expand vertically, displacing the older freshwater below (Fig. 4b). It appears to be a lagged response of the beach shortening event in February 2022 and the corresponding drop in groundwater head around ML3. The lower head propagates through the highly permeable glacio-fluvial deposits in landward direction and in turn increases the vertical hydraulic gradient over the tidal flat deposits, which, as a consequence, led to elevated seawater infiltration rates and expansion of the USP (Fig. 4a). The effect is also visible in the groundwater age of the seawater component (Fig. 5). At ML3–24 m, the freshwater water component was very old at the time, ranging from 18 to 25 years suggesting that the freshwater originated from the interior of the freshwater lens. At ML3–6 m, however, the freshwater component was younger, with an age ranging between 2 and 12 years, reflecting freshwater infiltration at the beach. With time, the groundwater ages increased around ML3, while correspondingly salinity decreased. This also appears to be a result of the freshening trend in response to the beach shortening event in February 2022. Groundwater age is/was most variable rather offshore, well seaward of ML3 (Fig. 6). Here, the oldest freshwater was up-coning and mixed dynamically with the very young seawater and younger freshwater components. Temporal variation in groundwater age was also visible at the outer fringe of the fresh-seawater mixing area in 20–30 m depths between ML1 and ML2 (Fig. 6). Here, young fresh and brackish water from the beach and old freshwater from the freshwater lens converge and thus increase the age gradient, so that a moving transition zone has the highest impact on groundwater age variability. Sediment accretion and erosion at the beach change the hydraulic gradient, also affecting greater depths at the upper beach. For instance, the February 2022 beach erosion event slowly pulled younger water, mainly freshwater, from the upper beach into greater depths. This is visible by the decrease of groundwater age from more than 17 years prior to the event to around 7 years afterwards (Fig. 5). Prior to the February 2022 event, the morphology of the beach was characterised by a ridge and runnel system which was subject to continuous migration up-and down the beach (see Video S1). This led to multiple in- and exfiltration locations that shifted over time. As consequence, strong lateral gradients in groundwater age and temperature developed in the shallow beach in shore-perpendicular direction, which were not position-stable. These gradients were visible even where salinity gradients were less pronounced, because of upcoming of older seawater with a damped temperature signal.

Groundwater temperatures were around 11 °C in large parts of the model domain (Fig. 4c). Temperature dynamics, mainly owing to seasonal changes, were confined to the shallowest part of the beach aquifer (Fig. 6). Temperature changes quickly dampened at greater depth, where the temperature stabilized (Fig. 3b and c). A curved, position-stable and continuous temperature band of around 8–9 °C stretched from the upper beach between ML1 and ML2 to below ML2–24 m (Fig. 4c). These lower local temperatures result from recurring winter storm floods that infiltrate into the upper beach and are rapidly transported downwards, promoted by the high hydraulic conductivities in this area (Fig. 2).

3.3 Fluxes and groundwater ages of the fresh and saline SGD components

The simulated fresh and saline SGD water fluxes vary on different time-scales and are in the range of 1.0–5 and 2.0–17 m² d⁻¹, respectively (Fig. 7a). High-frequency variations can be seen on daily time-scales, as well as on spring-neap time-scales. Variation on daily time-scales are attributed to meteoric effects altering the astronomical tidal range. Changes between more than 3 and 10 m² d⁻¹ for the fresh and saline component, respectively, can be observed (Fig. 7a). High-frequency variations in SGD water flux stemming from semi-diurnal tides variation and spring-neap tides were also found by Heiss and Michael (2014) for a beach in Cape Henlopen, Delaware, USA. The longer-term variations for the saline component (smoothed red line in Fig. 7a) appear to loosely correlate with the changes of spatial variation in beach morphology in terms of the beach slope's standard deviation, rather than the mean beach slope in the intertidal zone (Figs. 7c, S10a and b). Thus, the sudden drop in the saline SGD water flux at and after the beach shorting event in February 2022 seemed to be attributed to the abrupt decrease in beach slope variability rather than the step-increase in mean beach slope. The flux recovered as the beach slope variability recovered. Previous model-based saline SGD flux estimates from May 2014 (Beck et al., 2017) were about 4.3 m² d⁻¹. This is only a little bit lower than the here computed flux of about 5 m² d⁻¹, but otherwise agrees very well with the current model results. The lower flux computed by Beck et al. (2017) is expected, as their model considered fixed beach morphology, mean tidal range and did neither account for dynamic wave-set up nor storm floods. The variation of the fresh SGD flux is mainly attributed to the seasonality of groundwater recharge on the island.

The simulated flux weighted mean groundwater age of the fresh SGD component range from 4 to 25 years, with an average of around 15 years (Fig. 7b). Longer-term variation is also superimposed by high-frequency variations on semi-diurnal and spring-neap time-scales. Correlation with the spatial variation of the beach slope is not visible. The longer-term variability appears to depend to some degree on the saline SGD flux (Fig. 7a). It seems that predominantly younger freshwater originating from meteoric recharge at the upper beach is being pushed out during periods of elevated saline SGD fluxes. The flux weighted mean groundwater age of the saline SGD component varies between 0.75 and 2.3 years, and correlates with both, the spatial variation in beach morphology and the mean beach slope in the intertidal zone (Fig. S10c and d). The phenomenon of higher saline SGD flux and its younger mean age with increasing relief in the intertidal zone can be explained by so-called Tóthian flow. In his theoretically analysis, Tóth (1963) showed that topography-controlled flow with an increasing number of valleys and hills within a basin (and thus increased variability in the topographic gradient) leads to more local and shallower flow cells of higher discharge and shorter residence times.

Figure 7 (a) Simulated time series of saline (red) and fresh (blue) SGD water flux. Thin lines represent the fluxes derived half-daily after each stress period, thick lines represent the signal smoothed with a Savitzky-Golay filter with a polynomial order of 2 and a window length encompassing 101 data points. (b) Simulated time series of flux weighted groundwater ages of the saline (red) and fresh (blue) SGD components. Thin lines represent the fluxes derived half-daily after each stress period, thick lines represent the signal smoothed with a Savitzky-Golay filter with a polynomial order of 2 and a window length encompassing 101 data points. (c) Mean and standard deviation of the beach slope in the intertidal zone. Only the time-period from January 2013 onwards are presented, i.e., the period were the dynamics of the sea side boundary condition is considered in the model.

3.4 Limitations

In the presented field-site model approach, groundwater flow was assumed to be two-dimensional and lateral, longshore gradients were neglected. Other studies highlighted the importance of three-dimensional flow in beach aquifers (Zhang et al., 2016).

Further, simplifying assumptions on the morphology of the subtidal zone were made due to a lack of bathymetric data, although bathymetry is known to affect seawater intrusion rates and the interface of freshwater tube and saltwater wedge (Walther et al., 2017). Perhaps a more detailed description of the subtidal morphology and its changes in time could yield better calibration results at ML3–6 m and ML3–12 m. The assumed gradual change in intertidal beach morphology due to linear interpolation in time from 6-weekly profiles may not be adequate when sudden accretion or erosion events occurred. Nevertheless, at least during the present observation period, higher resolution profiles in time seemed unnecessary in order to describe the observed dynamics in the considered state variables.

The tide-averaged sea-side boundary condition is a drastic simplification of the complex dynamic seawater-groundwater exchange process in response to tides and waves. Calibration of a 20-year density-dependent groundwater flow and transport model on this spatial scale for a real field site, however, is virtually impossible when explicitly accounting for tides due to the large numerical burden. Interestingly, the observed changes in salinity and temperature, even in the shallower region of the beach aquifer, suggest that they are more controlled by average hydraulic gradients rather than by flow reversals and varying in- and exfiltration rates at the semi-diurnal timescale.

Variably saturated flow effects were neglected in the present model. Recent simulations of Luo et al. (2023) revealed that the effect of variably saturated flow decreases the phase averaged intertidal water flux by around 20 %. The authors, however, simulated an idealized, rather shallow aquifer of 15 m depth and thus a larger capillary-fringe-to-aquifer-thickness ratio as at the here studied beach aquifer, which promotes the importance of variably saturated flow. Even more recent, Luo et al. (2024) investigated the effect of variable saturated flow on redox zoning in a reactive transport modelling study with a similar model setup. The authors noticed some differences in the extent of the oxic zone when comparing simulation results with and without unsaturated zone effects. However, in the same study, benchmark simulations with two different codes, both considering unsaturated flow, yielded discrepancies of a magnitude similar to the effect of unsaturated flow. This fact and the generally known uncertainty of hydraulic conductivities, Genuchten-Mualem parameters, porosities, heterogeneity and geological layering, in combination with dynamic tidal-ranges, storm events, varying wave setup and morphodynamics at real field sites likely overrule effects that stem from variably saturated flow.

It needs to be pointed out that although rigorous model calibration was carried out, a perfect fit to observations was not expected nor intended. The aim was rather to obtain a model that is able to capture the principal patterns and trends observed in the field.

4 Conclusions

To our knowledge, this study is the first to develop a comprehensive density-dependent groundwater flow and transport model for a real high-energy beach aquifer that is calibrated to an extensive high resolution dataset. Despite the discussed model assumptions and simplifications, the model matched the observations very well. Model analyses revealed groundwater flow and transport patterns as well as information on subsurface residence times of the fresh- and seawater endmembers, and how the overall groundwater age was affected by mixing.

Further, hot spots of high temporal variability in salinity, groundwater age and temperature could be identified. The most dynamic fresh- and seawater endmember mixing was observed at the fringes of the infiltration and exfiltration zones of the USP. Here, mixing-controlled reactions (Meyer et al., 2025a), e.g., the oxidation of dissolved ferrous iron and precipitation of iron(hydr)oxides as so-called iron-curtain (Charette and Sholkovitz, 2002), or nitrification due to ammonium oxidation (Spiteri et al., 2008) may be promoted. Groundwater age dynamics are confined to a rather narrow belt below the USP, as well as to a zone further offshore, well seaward the low-tide mark. In these zones, biogeochemical turnover and secondary reactions may be controlled by the varying degradability of dissolved organic matter as a function of groundwater age (e.g., Waska et al., 2021). Temperature variations were only large in the uppermost part of the beach aquifer. Here, strong shifts and variable extents of redox-zones can be expected. Thus, the identified hotspots of temporal variability may also be reactive hots-spots. This means that our model can be used as a guide for field sampling campaigns targeting at subsurface biogeochemical transformation processes. Nevertheless, dense monitoring at high resolution and over a few years are needed to obtain a reasonably calibrated numerical model that allows to visualize a detailed picture of the flow and transport processes in the subsurface of a high-energy beach.

Further, the model results show that saline SGD fluxes vary considerable at daily and spring-neap time scales. The astronomical tide is always superimposed by meteoric effects which alter high water and low water levels from day to day. These are not easily predictable and therefore cannot be reflected in predictive modelling to guide field sampling campaigns on measuring these fluxes, e.g., by seepage meters. The model results suggest that a robust mean value can only be obtained when measurements capture these daily dynamics over a full spring-neap tide cycle. On longer-term, i.e., monthly to yearly time scales these fluxes gradually change and appear to be mainly controlled by the spatial variability of the beach slope rather than its mean gradient. The large dynamics in groundwater age of the fresh SGD component seem to loosely correlate with the magnitude of the saline SGD flux. Higher saline flux also pushes out younger fresh water that originates from meteoric recharge at the upper beach.

The calibrated model yields information on the permeability distribution within the beach aquifer. Whether the incorporation of new data and continuous re-calibration of the model confirms the present conceptual model remains to be investigated. The present model will serve as tool for further analyses. For example, the importance of individual boundary effects (e.g., dynamics in beach morphology, storm events, or wave-set) on mixing and residence times, which are critically important for understanding and quantifying reactive transport in the present subterranean estuary, will be investigated.

Appendix A

Table A1 Calibrated and fixed model parameters.