Abstract

The drinking water supply on barrier islands largely depends on freshwater lenses, which are also highly relevant for island ecosystems. The freshwater lens presented in this study is currently developing (since the 1970s) below the very young eastern part of the North Sea barrier island Spiekeroog, the so-called “Ostplate.” Due to the absence of coastal protection measures, formation, shape, and extent of the freshwater lens below the Ostplate are unaffected by human activities but exposed to dynamic changes, e.g., geomorphological variations and storm tides. The main aim of this paper was to reconstruct the evolution of the freshwater lens over several decades in order to explain the present-day groundwater salinity distribution. In addition, the study assessed the impact of geomorphological variations and storm tides on the freshwater lens formation. Detailed field observations were combined with a transient 2-D density-dependent modeling approach. Both field observations and simulations show an asymmetric freshwater lens after ~42 years of formation, whereby the horizontal extent is limited by the elevated dune area. The simulations indicate that the young freshwater lens has nearly reached quasi-steady-state conditions mainly due to the continuous mixing with seawater infiltrating during storm tides, which inhibits further growth of the freshwater lens on the narrow island. The findings further show that (i) a neglection of storm tides results in a significant overestimation of the freshwater lens extent, and (ii) the modeled present groundwater salinity distribution and shape of the freshwater lens are predominantly determined by the position and extent of the elevated dune area at the past ~20 years. Hence, annual storm tides have to be directly implemented into numerical models to explain the groundwater salinity distribution and the extent of young freshwater lenses located in highly dynamic tidal environments.

1. Introduction

The water supply in coastal regions largely depends on groundwater from local aquifers. Since almost 25% of the global population lives within 100 km of a coastline, the high demand for freshwater creates an immense pressure on the coastal groundwater resources [13]. With the expected doubling of the population density in coastal regions by 2030 [4], the demand for freshwater and its extraction will further increase. The available freshwater resources in coastal areas and related ecosystems are highly vulnerable and often threatened by salinization due to seawater infiltration and overuse. The vulnerability of coastal aquifers is additionally enhanced by saltwater intrusion and more frequent inundations due to increasing storm frequencies both coming along with climate change and predicted sea level rise [1, 5, 6].

Barrier islands are commonly characterized by underlying fresh groundwater resources, which typically exist in the form of freshwater lenses (FWLs) [7, 8]. On dune islands, like the barrier islands located in the southern North Sea, the formation of dunes allows for the development of FWLs due to their elevation and infiltration capability. Reaching heights above the maximum level of storm tides, elevated areas are no longer prone to inundations, and the previously saline unconfined aquifer gets recharged by precipitation, displacing the saltwater and reducing the groundwater salinity [911]. Due to density differences between fresh and saline groundwater, the freshwater floats on top of the underlying saltwater (e.g., [1216]). The transition zone that separates fresh and saline groundwater normally is not a sharp boundary as assumed by analytical solutions [7, 8]. It is rather a zone of up to several meter thickness due to mixing and hydrodynamic dispersion, displaying an increase in salinity with depth (e.g., [10, 13, 17, 18]).

Both the FWL and the transition zone are dynamic and influenced by natural factors including tides, seasonal changes in groundwater recharge/discharge, ground surface elevation, and island stratigraphy [14, 1721]. In the absence of coastal measures like dykes, inundation events can cause salinization of the FWL, posing a serious threat to local water supply [2226]. The inundation frequency is, thereby, mainly determined by ground surface elevation above sea level [27]. Given the fact that sea level and the frequency of storm tides are predicted to rise, barrier islands and the available freshwater resources, respectively, will become even more vulnerable to inundations and saltwater intrusion, especially on islands where topography constrains the rise of the FWL [5, 28, 29].

According to Illangasekare et al. [30], the infiltration of saltwater through the unsaturated zone during inundation events has the largest spatial impact on the FWL. Its impact on the FWL depends, thereby, on the thickness of the unsaturated zone, i.e., the thicker the unsaturated zone, the larger the amount of saltwater intruded [30, 31]. The FWL recovery after salinization due to storm tides to preinundation levels can range between months [2, 22, 32] and several years [23, 25, 33]. It is mainly driven by recharge from precipitation flushing out the saltwater [30], and therefore, the time of recovery is a function of the recharge rate [34, 35]. Besides the thickness of the unsaturated zone, hydraulic properties of the aquifer, i.e., hydraulic conductivity, and geological heterogeneities also exert control on the recovery of the FWL [22]. For instance, Yang et al. [36] conducted a modeling study and showed that a low horizontal permeability and a high vertical permeability increase the aquifer vulnerability during inundation events. Infiltrating seawater may appear as descending saltwater fingers due to density differences between the fresh- and saltwater, causing a rapid downward migration of the intruded saltwater after inundation events resulting in an acceleration of the recovery of the upper parts of the FWL, but a long-lasting deterioration of the deeper parts [22, 25, 32, 33, 3739].

Furthermore, overextraction of fresh groundwater resulting in saltwater intrusion but also contamination with biological and chemical pollutants from surface sources are anthropogenic factors affecting thickness, quality, and quantity of FWLs [18, 21, 24, 40].

A FWL is currently developing below the eastern part of the North Sea barrier island Spiekeroog, the so-called “Ostplate,” representing the study area. The formation of the FWL started in the 1970s, when the first dunes covered with Ammophila arenaria (European beachgrass) appeared [41]. As a nature conservation area located within the national park “Nationalpark Niedersächsisches Wattenmeer” and thus part of the World Heritage Site “Wadden Sea,” declared in 2009, the Ostplate remains largely unaffected by human activities. The island’s village and infrastructure, including the groundwater abstraction wells, are all located in the western part of Spiekeroog, whereas the younger Ostplate and the formation of the associated FWL are unaffected by anthropogenic activities. Coastal protection measures do not exist, and the FWL is exposed to natural dynamic changes like storm tides and geomorphological variations. The FWL of the Ostplate has so far not been monitored, and except for preliminary estimates from numerical simulations by Röper et al. [41], vertical extent, tidal and seasonal dynamics, and temporal evolution of the FWL are yet unknown. However, Holt et al. [27] reported that the shallow groundwater salinities (groundwater salinity near the water table) are influenced by storm tides and inundation frequency. The near-surface extent of the FWL of the Ostplate is therefore characterized by spatial and temporal variations.

Processes affecting FWLs were analyzed in many numerical modeling studies. These studies cover natural and anthropogenic influences regarding the (formation of) FWLs (e.g., [18, 21]) and response and recovery of FWLs after saltwater intrusion due to inundation with seawater (e.g., [25, 26, 32]) or are generic model approaches investigating impacts of potential and future storm tide and inundation events (e.g., [22, 31, 42]). Thereby, a range of simulations conducted on islands are restricted in the sense that they only take one sea boundary into consideration, i.e., the model domain does not cover the entire island width, that effects of storm tides on the formation of FWLs are neglected, or only single storm tide and inundation events are considered. In contrast, Huizer et al. [43] showed for an artificial environment (beach nourishment) at the coast of the Dutch mainland how storm tides and geomorphological changes affect the formation of a very young FWL.

In the present study, we combine detailed field observations of a currently developing FWL at a barrier island with a comprehensive modeling approach. In contrast to a previous study at the site [27] that focused on the spatial and temporal dynamics of shallow groundwater salinities as a function of inundation frequency, the aims of this study were (i) to reconstruct the evolution of a currently developing FWL impacted by tides and storm tides over several decades in order to explain the present-day groundwater salinity distribution within the aquifer and (ii) to assess the impact of geomorphological changes and storm tides on the FWL formation.

2. Study Area

The barrier island Spiekeroog is part of the East Frisian Island chain, located ~6.5 km in front of the northwest coastline of Germany in the North Sea. It lies between the islands of Langeoog to the west and Wangerooge to the east. The west-east and north-south extents of the present island are ~9.8 km and 2.3 km, respectively, resulting in a total area of the island of ~21.3 km2. Whereas the older western part of the island (~4.5 km) is inhabited, the younger eastern part (~5.3 km), the Ostplate (Figure 1), is uninhabited and highly protected by nature conservation regulations [44, 45].

2.1. Geomorphology of the Ostplate

Around 1940, the Ostplate was no more than a periodically flooded sand flat. The development of the Ostplate towards a dune island was caused by land reclamation and embankment activities at the mainland [41]. According to Röper et al. [41], the dune formation started with the initial growth of unstable pioneer dunes due to aeolian processes in the 1940s. Using Optically Stimulated Luminesce (OSL) age dating on core sediments, which were extracted at the dune base, Seibert et al. [46] could confirm that the aeolian sediments of the upper dune parts are younger than ~80 years. The first appearance of A. arenaria around 1970 marked the formation of white dunes (secondary dunes) and the start of the freshening process in the dunes [44, 47]. With a maximum height of ~11 m above sea level (asl), the present dunes of the Ostplate are mainly characterized as white dunes and partly grey dunes (tertiary dunes) [41]. The shelter of the dunes also enabled an extensive growth and spreading of salt marsh vegetation [27, 41]. For more detailed information on the evolution of the Ostplate, see Holt et al. [27], Röper et al. [41], and Seibert et al. [46].

2.2. Hydrology and Hydrogeology

Spiekeroog is affected by semidiurnal tides, with the tidal regime (mean tidal range of 2.72 m) being mesotidal according to the classification of Hayes [48]. Mean High Water (MHW) reaches 1.39 m asl and Mean Low Water (MLW) 1.33 m below sea level (bsl; annual means from 1995 to 2015 [49]). The mean spring range is 3.09 m, whereas the mean neap range is 2.23 m [50]. Three storm tides, reaching a maximum height of 3.05 m asl, were recorded during data collection in the winter 2016/2017 [49].

Affected by a temperate climate with all-season rainfall, the average annual precipitation amounts to 808 mm at Spiekeroog (from 1984 to 2011 [51]). With respect to the historical island formation [52], the evolution of the FWL in the western part of Spiekeroog started ~350 years ago, and the present FWL has a thickness of ~44 m and is being used for the island’s drinking water production [53]. Based on apparent 3H-3He ages, Röper et al. [53] calculated groundwater recharge rates of 300-400 mm a-1 for the western dune area of Spiekeroog, while other references reported on recharge rates between 301 and 640 mm a-1 [54, 55]. The salt marsh, in contrast, is characterized by lower recharge rates due to the less permeable fine-grained sediments (clay, silt) and drainage due to the existence of tidal creeks (as suggested by [56, 57]).

The sediments of the western part of Spiekeroog consist of fluviatile Pliocene sand deposits partly with enclosed peat layers, which are overlain by sandy glaciofluvial Pleistocene deposits and fine- to coarse-grained Holocene sands. At a depth of 40 to 60 m bsl, the sand deposits are underlain by a clay layer of 1.5-15 m thickness [55] (compare [53, 58]).

Röper et al. [41] also proved the existence of freshwater below the dunes at the Ostplate, and Holt et al. [27] could show that several separated FWLs, restricted to the elevated dune area, exist here (Figure 1). Strongly affected by spatial and temporal water-level variations, most importantly the storm tides mainly occurring in winter, the largest near-surface extent of the FWLs was observed in summer, whereas the winter storm tides lead to a shrinkage. The northern beach area and the broad salt marsh covering most of the southern area of the Ostplate are characterized by underlying brackish to saline groundwater as a consequence of frequent inundations [27]. However, a reconstruction of past shallow groundwater salinities indicates an expansion of areas underlain by freshwater due to dune growth and a simultaneous saltwater decrease for the past 15 years [27]. Seibert et al. [46] could further show that the groundwater salinity at the margin of the here investigated FWL also shows high temporal variations in greater depth and is largely determined by storm tides. Seawater infiltration during inundations and salinization of the former fresh margin of the FWL results in an abruptly steep rise of the salinity followed by a gradual salinity decrease over time.

3. Material and Methods

3.1. Instrumentation and Sampling

Six monitoring sites were installed in a north-south trending transect at the western part of the Ostplate (A-A’, Figure 1) in March 2016, to facilitate groundwater level and electrical conductivity (EC) measurements. The transect is oriented along the principle assumed north-south groundwater flow direction that follows the general topographic gradient (Figure 1). Multilevel observation wells (OWs) were installed at the beach (monitoring site B, 3 piezometers with a single screen each; B-a: 0.83 m asl, B-b: -5.31 m asl, and B-c: -12.69 m asl) and at the northern (monitoring site DN, 4 piezometers with a single screen each; DN-a: 1.54 m asl, DN-b: -1.54 m asl, DN-c: -4.55 m asl, and DN-d: -6.44 m asl) and southern (monitoring site DS, 3 piezometers with a single screen each; DS-a: 1.70 m asl, DS-b: -1.04 m asl, and DS-c: -3.07 m asl) dune base, respectively (Figure 2). Three shallow single screen OWs were installed within the salt marsh south of the dunes (Figure 1), screening near-surface groundwater at 1.29 m asl (SM-1), 1.34 m asl (SM-2), and 0.97 m asl (SM-3) (Figure 2). All of the installed piezometers have a screen length of 1 m.

Groundwater levels were continuously measured at all OWs (except for B-b and DN-b) with CTD-Divers (DI271, Schlumberger Water Services) and with Mini-Divers (DI501, Schlumberger Water Services), respectively, using 10-minute intervals from April 2016 to August 2017. The EC was measured manually during six hydrochemical sampling campaigns (bimonthly) described in Seibert et al. [46]. Groundwater samples were extracted from all OWs (Figure 2) using a submersible pump (Gigant, Eijkelkamp). In addition, the Direct-Push method (probe: water sampler DP32/42 with 1 inch internal filter of 0.5 m length, Stitz) was applied at the monitoring sites BN, DN, and DS in September 2016, which allowed for the depth-specific extraction of additional groundwater samples () and a higher spatial resolution of the vertical FWL extent compared to the permanent monitoring sites DN and DS alone. Samples were collected every meter (Figure 2) using a foot valve pump (S-60, Stitz). Table S.1.1 (see Supplementary Material) gives an overview of the instrumentation of the installed OWs and Direct-Push drillings as well as time periods of the water level and EC measurements.

The EC of all samples was measured using a Hach HQ40d multidevice (automatic temperature correction to 25°C) and subsequently converted into salinity using the method described in Holt et al. [27]. Groundwater salinities were classified according to Freeze and Cherry [60], with freshwater, brackish water, and saltwater having a salinity of <1, 1-10, and >10 (g L-1), respectively.

3.2. Numerical Model

Assuming the principle north-south groundwater flow direction (see Section 3.1 and Supplementary Material S.2), a transient vertical density-dependent 2-D groundwater flow model was set up for the cross-section A-A’, i.e., along the transect of OWs (Figure 1), using the software SEAWAT [61].

The north-south oriented model cross-section A-A’ (Figures 1 and 3) is bounded at the MHW-mark (1.39 m asl) north and south of the Ostplate and includes the beach in the north (~520 m length) with a slope of about 0.002, as well as dune (~100 m length) and salt marsh areas in the south (~1,610 m length). The extent of the model is 2,230 m in the horizontal direction, with a discretization of 10 m. The topography was obtained from a laser-scan elevation model of 2014 (provided by the NLWKN [59]), interpolated with the field interpolator of the groundwater simulation software PMWIN (Processing Modflow, Version 8.0.47, www.simcore.com) to the model grid, and assigned to the top of the uppermost model layer. Except for the first two model layers, the vertical discretization is 1 m. The thickness of layer 1 ranges from 0.75 m to ~8 m, and the thickness of layer 2 ranges from 0 m to 0.9 m. This was done to circumvent the SEAWAT-specific rewetting process of dry cells during groundwater table movement across model layers, which sometimes causes numerical problems.

A no-flow boundary was assigned to the bottom boundary, assuming that the clay layer detected at the western part of the island (see Section 2.2) also exists at the eastern part of the island. Because of the extremely low morphological gradients of the beach and salt marsh areas and the correspondingly extremely small time steps (in the order of seconds) that would be needed to adequately resolve the extremely fast pressure propagation of the tidal signal at the top boundary, tide-averaged constant heads rather than real tides (following the approach of Vandenbohede and Lebbe [11, 62]) were defined at the northern and southern vertical model boundaries at the MHW-mark (Figure 3). The assigned heads were based on field measurements near the MHW-line at the beach at the western part of Spiekeroog [63] and slightly calibrated to 0.8 m asl (northern boundary) and 0.9 m asl (southern boundary). Reflecting the so-called tidal overheight due to asymmetric in- and outflow conditions during high and low tide, respectively, at a sloping beach [64], the hydraulic heads at the MHW-line are elevated compared to the mean sea level. A direct implementation of the tides (e.g., semidiurnal, spring, neap, and storm tides) would lead to excessively large computing times (similar to, e.g., Comte et al. [65]) and is to the author’s knowledge not feasible for a simulation of >40 years. Furthermore, a detailed analysis of the daily tides would be beyond the scope of this study, as it is focusing on the FWL evolution on a decadal time scale. The northern and southern vertical boundaries for the transport model were defined as non-dispersive flux boundaries [66]. Here, the salinity of the water entering the model domain through the constant head flow boundary was set to seawater salinity (32 g L-1), whereas for outflow conditions, it was set to the salinity of the water computed at the boundary cell. The non-dispersive flux boundary is commonly applied to simulate seawater-groundwater interaction at the costal boundary (e.g., [67, 68]).

The horizontal hydraulic conductivity served as a calibration parameter, and the final calibrated value of agreed well with the hydraulic conductivity derived from grain size analysis and Darcy experiments (see Supplementary Material S.3, Fig. S.3.1). It is within the range of values for sand [69] and in line with results of pump tests conducted at the western part of Spiekeroog ( and ) [55]. The vertical anisotropy of the hydraulic conductivity was set to 1 : 1, as anisotropic hydraulic conductivities yielded poorer model results (similar to Schneider and Kruse [21]). The effective porosity was set to 0.35, which is a similar value found in related studies (e.g., [10, 11, 41, 62]). The longitudinal dispersivity and the vertical transverse dispersivity were and , respectively, and found during model calibration. The model parameters are listed in Table 1.

Assuming that the FWL formation started around 1975 (based on Röper et al. [41]), the total simulation time was 42 years (1975-2016).

The conceptual model described above was used to explain the present-day groundwater salinity distribution and extent of the FWL, and to reconstruct a temporal evolution of the FWL below the dune belt in the western part of the Ostplate. Due to missing historical data, the simulation is only calibrated with respect to the present groundwater levels and groundwater salinities (2016-2017). The simulation assumed at first a constant geomorphology, i.e., only the geomorphology of the year 2014 was considered. This base case is referred to as Constant Geomorphology-model (CG-model) in the following. Furthermore, two additional model cases based on the CG-model were considered. The first additional case indirectly accounted for geomorphological changes (in the following referred to as Variable Geomorphology-model (VG-model)). Within the second additional case, the CG-model was run without implementing annual storm tides (in the following referred to as Constant Geomorphology Without Storm Tides-model (CGWST-model)). Comparing both additional cases to the base case (CG-model) allowed to study how the FWL and the groundwater salinity distribution depend on geomorphological variations and storm tides that are characteristic for the highly dynamic environment investigated here. We would like to emphasize here that the simulations are not intended to exactly replicate the observed hydraulic heads and groundwater salinities over time at the cross-section investigated here which is not possible for the lack of data before 2016. Rather, the simulations were conducted to reflect the observed pattern of hydraulic heads and groundwater salinities over time in a more generic, however, real-data oriented sense, in order to analyze the significance of a varying geomorphology and storm tides on the evolution of a FWL on a young barrier island.

The groundwater flow was numerically solved with the Geometric Multigrid Solver (GMG) and advection-dispersion with the Hybrid Method of Characteristics (HMOC) in conjunction with the Generalized Conjugate Gradient (GCG). The HMOC scheme provides a solution with virtually no numerical dispersion independent from the Grid Peclet number [66] but is prone to small mass balance errors. In all our simulations, a mass balance error of <1.5% as well as grid-independent of the simulation results were ensured.

3.2.1. Constant Geomorphology with Storm Tides: CG-Model

A specified flux boundary for the flow model with different recharge rates was assigned along the top boundary, coupled with a non-dispersive flux boundary for the transport model characterized by the salinities of the recharge water (shown by the salinity distribution of annual recharge in the last row from above labelled “2014” and by the salinity distribution of recharge during storm tide in Figure 3). Recharge rates were considered as calibration parameters. Nevertheless, the calibrated rates correspond to the sediment type along the cross-section A-A’, the unsaturated zone thickness and the geomorphology of the year 2014. Due to the lack of historical data prior to 2014, rates were applied to the entire simulation time (1975-2016). The highest recharge rates were assigned to the dune area (fine to coarse sand, recharge zone C, Figure 3 and Table 2) and a 30 m wide area north of the dunes. The beach was divided into two recharge zones (recharge zone A and B, Figure 3 and Table 2). Although other studies conducted on Frisian Islands reported on lower recharge rates (e.g., [25, 56, 70, 71]), the range of recharge rates used in this study is in line with previously determined rates at Spiekeroog (see Section 2.2). Salt marsh areas were divided into three recharge zones (recharge zone D, E, and F, Figure 3 and Table 2) with lower recharge rates as those areas are characterized by less permeable near-surface fine-grained sediment layers (i.e., clay and silt), which were also described in Sulzbacher et al. [56], as well as by draining tidal creeks. Also, a denser vegetation cover occurring in salt marshes typically diminishes recharge due to enhanced transpiration [21].

Fresh recharge (salinity 0 g L-1) was implemented in the elevated dune area (≥3 m asl). The recharge salinity increases with increasing distance from the dunes and decreasing ground surface elevation at the beach and salt marsh areas (shown by the salinity distribution of annual recharge in the last row from above labelled “2014” in Figure 3). The implementation of this salinity gradient enabled an indirect but efficient way of accounting for an increase of inundation frequency and its bulk effect on introduced salt mass with increasing distance from the dunes. Again due to the lack of historical salinity measurements, we applied the present-day near-surface groundwater salinity distribution [27] for the entire simulation period.

Severe flooding events inundate the entire Ostplate except for the elevated dune areas [27] roughly two times per year on average [49] (tides reaching a height >1.40 m above MHW, referred to as storm tides in the following) affecting the groundwater salinity distribution at the margin of the FWL [46]. To account for these, transient recharge rates (instead of recharge rates, we refer to infiltration rates for seawater infiltration in the following) and recharge salinities were implemented into the model. Volumes of infiltrating seawater were estimated based on the thickness of the unsaturated zone and the porosity (0.35), presuming that the entire unsaturated zone of inundated areas gets filled up with seawater during an inundation event (as suggested by Illangasekare et al. [30] and Chui and Terry [31]). The thickness of the unsaturated zone was calculated for the locations of the monitoring sites and linearly interpolated in between. Consequently, the highest infiltration rates were assigned to the northern and southern dune base.

Based on the two storm tides per year on average, we assumed that the unsaturated zone gets filled up with seawater twice a year. For numerical reasons, i.e., a long computing time due to small computational time steps and the vulnerability of the simulation to numerical instabilities, we distributed the respective infiltration volume evenly over a period of six days. Since salinities directly at the margin of the FWL were, however, overestimated (results not shown), seawater infiltration rates were subsequently reduced to one third in areas characterized by ground surface elevations >2.40 m asl. This is a reasonable assumption as these areas are submerged in seawater only during very short periods of time. Contrary to Chui and Terry [31], we assumed that this time is insufficient to fill up the entire unsaturated zone. The non-inundated part of the dune area is characterized by an infiltration rate of 0.52 m a-1 at all model times. Because of low permeabilities and a denser vegetation cover, only half of the calculated infiltration volumes were assigned at the salt marsh areas during storm tides.

Recharge salinities during storm tides correspond to local seawater salinity (32 g L-1), except for the non-flooded parts of the dune area (0 g L-1, shown by the salinity distribution of recharge during storm tide in Figure 3).

Table 2 gives an overview of the implemented infiltration/recharge rates, the total infiltration/recharge, and the corresponding recharge salinity distribution during storm tides as well as during the remainder of the year, respectively.

3.2.2. Variable Geomorphology with Storm Tides: VG-Model

Temporal changes of the geomorphology at the Ostplate were evident based on (a) the analysis of aerial and infrared pictures, and digital elevation models for the years 1985, 1991, 1999, and 2011 by Röper et al. [41] and (b) digital elevation models of the years 2008 and 2014 (provided by the NLWKN [59]). We presume that these geomorphological variations influenced the formation of the FWL. Areas with fresh recharge relevant for the FWL formation (dune area ≥3 m asl) were located for the years 1985, 1991, 1999, 2008, 2011, and 2014 using ArcGIS (see Supplementary Material Fig. S.5.1).

In cases where the time intervals between the mentioned years is >5 years, the extent of the dune area was linearly interpolated (Table S.5.1), resulting in interpolated extents for the years 1988, 1995, and 2003. The extent of the dune area was set to 10 m, as the minimal extent of a model cell, for the year 1975 (assumed start of the FWL formation [41]) and to the half of the dune extent of 1985 for the year 1980.

The implementation of the temporally varying geomorphology into the model set-up was done by adjusting the recharge zones and salinities, respectively, which were identified and used in the CG-model, while the remainder of the model parameters (Table 1) was left unchanged. Figure 3 gives a schematic overview of the changing geomorphology implemented indirectly in the model by an altering domain of fresh recharge (blue color) for the aforementioned years (shown by the salinity distribution of annual recharge). For simplicity, the storm tide salinity distribution is displayed for the year 2014 only, but also changed depending on the geomorphology for the years before (Figure 3). The same applies to the location and extent of the recharge zones. Recharge and infiltration rates, respectively, and recharge salinities depending on the position and extent of the dune area in the respective years are listed in Table S.5.2 and Table S.5.3 (see Supplementary Material).

3.2.3. Constant Geomorphology without Storm Tides: CGWST-Model

Apart from the infiltration rates and recharge salinities assigned for the storm tides in the CG-model, the model parameters were the same as in the CG-model (Tables 1 and 2).

4. Results and Discussion

4.1. Present-Day Groundwater Salinity Distribution and Extent of the FWL: Field Data

The schematic cross-section shown in Figure 2 gives an overview on the groundwater salinity distribution along A-A’ and shows the results of the sampling campaign in September 2016. Freshwater (salinity ≤1 g L-1) was only found in the upper OWs and Direct-Push locations below the dunes (DN-a, DN-3,4, DN-b, DS-a, and DS-4,5), indicating a FWL with an asymmetric shape. The change from fresh (salinity ≤1 g L-1) to brackish groundwater (salinity 1-10 g L-1) occurs at ~5 m below ground surface (m bgs) at both the northern (between DN-b and DN-6) and southern dune base (between DS-5 and DS-b). Hence, the maximum vertical extent of the FWL is ~4 m at the northern margin of the FWL (related to the mean depth between DN-b and DN-6). Compared to the FWL thickness of 14-25 m bsl suggested by preliminary numerical simulations of Röper et al. [41], the actual FWL extent is significantly smaller.

The brackish transition zone below the FWL is not a sharp boundary but a diffuse zone of several meter thickness with increasing salinities with increasing depth (DN-6, DN-c, DN-d, DS-b, and DS-c). However, the exact vertical extension of the brackish zone remains unclear, since it was not possible to drill deeper and saline groundwater was not intersected even at the deepest wells at the dunes (i.e., DN-c, DN-d).

The storm tides in December 2016 and January 2017 resulted in inundations of the Ostplate leading to a shrinkage of the FWL at the northern margin, and freshwater was only found at the uppermost OW south of the dunes in January 2017 (DS-a, see Supplementary Material Fig. S.6.1).

4.2. Model Results

The groundwater salinity distribution in different years after the onset of the freshening process is illustrated in Figure 4 for all three model cases. Videos of the evolution of the FWL are available in the Supplementary Material (videos S.7.1, S.7.2, and S.7.3).

4.2.1. Temporal FWL Evolution until Present with and without Geomorphological Variations

Both the CG- and the VG-model show the development of an asymmetric FWL which horizontal extent is limited by the elevated dune area (≥3 m asl), i.e., the area of fresh recharge (CG- and VG-model in Figure 4). The asymmetrical shape of the FWL is caused by the location of the water divide south of the dunes in the upper part of the salt marsh at a distance of about 1,000 m (CG-model) and 980 m (VG-model) from the northern model boundary, respectively. Hence, all freshwater that infiltrates in the dunes flows towards the beach and the northern MHW-boundary, respectively. Initial freshwater bodies (salinity ≤1 g L-1) developed after 11 and 25 years in the CG-model and the VG-model, respectively. The time lag between the two model approaches is due to the strong movement of the dunes in the first simulation years, which was only accounted for in the VG-model.

The simulated maximum vertical extent of the present FWL is about 3-4 m, which is well in line with the vertical thickness observed in the field (Figure 2). The horizontal extent varies between 50 m (CG- model) and 40 m (VG-model). Considering a west-east extent of 1 m and a porosity of 0.35, the volume of the FWL is estimated to be ~65 m3 (CG- model) and ~38 m3 (VG-model), respectively, along the cross-section A-A’. Assuming one connected FWL below the dunes on the entire Ostplate with a west-east extent of ~3,600 m and a constant width of the dune ridge, the freshwater volume would amount to ~234,000 m3 (CG-model) and ~137,000 m3 (VG-model). The dimensions of the freshwater volume calculated here deviates strongly from the preliminary numerical simulations of Röper et al. [41], who calculated a freshwater volume ranging between 42 and 92 million m3. Reasons for the larger freshwater volume of their model may be the combined effect of assumed lower horizontal and vertical hydraulic conductivities, a wider area of fresh recharge, and the omission of storm tides.

The formation of the FWL is accompanied by the establishment of a wide brackish transition zone. Similar to field observations, the transition zone is not a sharp boundary but a diffuse zone of >10 m thickness showing increasing salinities with increasing depth (CG- and VG-model in Figure 4). The transition between freshwater and brackish water is located at ~5 m bgs (CG-model) and ~4 m bgs (VG-model). The distinctly thicker transition zone, compared to the vertical extent of the FWL, is in line with White and Falkland [18], after which FWLs characterized by a vertical extent of <5 m often are underlain by a transition zone larger than the FWL itself. The thickness of the transition zone is influenced by dispersive mixing, tides, temporal variabilities of recharge, hydraulic properties of the aquifer, and inundations evoked by storm tides [17, 18, 72].

Both simulation outcomes clearly suggest the appearance of saltwater fingering (onset in the CG-model: after 19 years, VG-model: after 26 years), especially below the beach due to infiltrating seawater. The development of saltwater fingers migrating from the upper beach towards the northern boundary at the MHW-line is caused by density differences between the groundwater characterized by lower salinities than seawater and the seawater which infiltrates in particular during storm tides. Such a behavior in the intertidal zone has previously been proposed by modeling and laboratory studies [39, 7375]. Saltwater fingers are also recognizable at the northern margin of the FWL.

To assess the model performance, modelled groundwater salinities and hydraulic heads (point water heads) after 42 years were compared to field observations from 2016, as no historical exists. The calculated heads of both models agree well with the observations, and the Nash-Sutcliffe model efficiency coefficient and the Root Mean Square Error (RMSE) only differ slightly between both models (Nash-Sutcliff: 0.93-0.94; RMSE: 0.03-0.04 m, see Supplementary Material Fig. S.8.1).

The comparison of mean calculated (last five years of simulations, see Supplementary Material Fig. S.9.1) and mean observed salinities (bimonthly sampling campaigns, see Section 3.1 and S.1.1) also shows a reasonable correlation for both models, despite some over- and underestimations. Again, both models exhibit only minor differences regarding the model fit with Nash-Sutcliffe coefficients of 0.90 and 0.91 for the CG-model and the VG-model, respectively (Figure 5). Since the saltwater fingering flow causes the groundwater salinity to vary substantially over time on the beach, a direct comparison between calculated and observed salinities for the year of sampling was not appropriate (similar to Post and Houben [25]). Particularly, the upper OWs at the northern margin of the FWL (DN) are characterized by a wide range of observed salinities (compare horizontal error bars in Figure 5). This is due to the highly dynamic environment, most importantly storm tides during winter and the associated salinization. Similar findings were also reported by Huizer et al. [76] for a beach nourishment site. The steep rise of the groundwater salinity following storm tides and the subsequent freshening which was also reported by Holding and Allen [22] for recent sedimentary islands, could be simulated (see Supplementary Material Fig. S.9.1, CG- and VG-model: DN-a) and are also visible in the calculated salinities at these OWs (vertical error bars in Figure 5, CG- and VG-model).

Apart from rather small differences in the spatial extent of the present FWL, the comparison of both models shows that the overall pattern of the current groundwater salinity distribution is very similar (CG- and VG-model in Figure 4). Modeled present-day groundwater salinity distribution and shape of the FWL are predominantly determined by the position and extent of the elevated dune area during the past ~20 years, when only slight variations of the fresh recharge area occurred (Figure 3). Geomorphological changes at the early stage of the FWL formation (VG-model, Figure 4) have only little impact on the present-day situation. Since a consideration of the geomorphological changes did not result in a better model fit, the simulation outcomes suggest that the geomorphological variations are not crucial for the present-day situation of the here investigated FWL. Hence, an overestimation of the FWL volume in the CG-model neglecting geomorphological changes, as reported by Huizer et al. [43], is not expected. This is because strong geomorphological changes only appeared at the beginning of the FWL formation, and the fresh recharge area has since continuously grown, while in the study of Huizer et al. [43], erosion was reported on.

4.2.2. Effect of Storm Tides on Groundwater Salinity Distribution and FWL Development

In contrast to the temporal variability of the geomorphology, the influence of storm tides is very pronounced. The comparison between the cases with (CG-model) and without (CGWST-model) the effect of storm tides clearly indicates that FWLs are very sensitive to storm tides and associated inundations (Figure 4). Furthermore, it shows that the lateral and vertical extent of the FWL largely depends on saltwater infiltration during inundation events. Similar to the findings of, e.g., Anderson Jr [23], Holding and Allen [22, 35], Huizer et al. [76], and Post and Houben [25], the vertical infiltration of seawater caused by storm tides resulted in severe groundwater salinization. It was not possible to achieve a satisfying model fit between observed and modelled data without accounting for storm tides. When neglected, simulated groundwater salinities were considerably underestimated, which especially holds for the salinities at the margin of the FWL and within the transition zone (Figure 5, CGWST-model). Moreover, the neglection of storm tides resulted in a significant overestimation of the FWL extent 42 years after formation (Figure 4), with a maximum vertical and horizontal extent of ~13 m and 180 m, respectively. Thus, considering storm tides is crucial for an adequate simulation of a long-term evolution of a FWL on a whole island width. Otherwise, the amount of infiltrating seawater impacting the development of a FWL would be underestimated, which was also reported by Huizer et al. [43] but neglected by Röper et al. [41]. However, the situation is likely to be different on older islands with higher morphology not prone to inundations.

Similar to the CG-model, saltwater fingers occurred that moved from the upper to the lower beach. Respective salinities were lower due to mixing with freshwater flowing from the larger FWL towards the sea and the missing salt mass from seawater infiltrating during storm tides (Figure 4). In the CGWST-model, the appearance of the fingering was caused by the tidal beach inundation inducing seawater infiltration which was simulated by the increasing recharge salinity towards the northern MHW-line.

4.2.3. Time of Equilibrium and Future FWL Development

The FWL in the CG-model approaches nearly steady-state conditions after ~35 years of formation (with regard to cells with salinity ≤1 g L-1; video S.7.1), while the FWL of the VG-model reaches nearly steady-state conditions after ~70 years of formation (results not shown). Further growth of the FWL is inhibited by the repeated infiltration of seawater caused by storm tides, and its final extent is very similar in both models after steady-state conditions are reached. Note that steady-state conditions in the VG-model would only be realized when the geomorphology is stagnant in the future. The extent of the ≤50% seawater area reaches steady-state after ~100 years of formation (CG-model, Figure 6), i.e., somewhat later than the FWL. Even the larger zone characterized by salinities ≤20 g L-1 does not expand further after 100 years due to downward migrating saltwater fingers.

However, transient conditions of the geomorphology, e.g., further dune growth due to accretion or dune erosion—which are likely to occur in such a highly dynamic environment studied here—will most likely affect the extent of the investigated FWL in the future, since variations of the coastal landscape have impact on the shape and volume of a FWL [77, 78]. This was also visible in the onset of the FWL formation in the VG-model (Figure 4). It has been shown previously that the topography of coastal landscapes significantly influences groundwater salinization of coastal aquifers following storm tides [79]. Linked to climate change and expected future sea level rise in the southern North Sea [5], length and frequency of storm tides and resulting inundations with seawater as well as the amount of fresh recharge are further constraints on the future FWL development [56, 80]. Moreover, fresh recharge depends on the vegetation cover in the dunes. As the dunes develop, vegetation will presumably increase and become denser, thereby causing a decrease in the recharge rate, which would consequently affect the FWL development [19, 65, 81]. Hence, real steady-state conditions will likely never be reached.

5. Conclusions

This study provides a comprehensive overview of the evolution of a FWL below a currently developing young barrier island located in a highly dynamic mesotidal environment unaffected by anthropogenic activities. Numerical simulations allowed to reconstruct the temporal evolution of the FWL and to assess the effect of storm tides and geomorphological variabilities on FWL formation. The following conclusions can be made: (1)Barrier island formation is accompanied by FWL formation as soon as the growth of dunes exceeds storm tide levels(2)Fingering flow appears to be a likely feature below beaches subject to freshwater discharge from FWLs in combination with frequent seawater inundations(3)For highly dynamic tidal environments subject to annual storm tides, these have to be directly implemented into models in order to explain the groundwater salinity distributions and the extents of young FWLs(4)Accounting for flooding frequency indirectly by adjusting the salinity of the recharge water accordingly appears to be a suitable method to mimic highly transient tidal systems without explicitly incorporating tides, as shallow groundwater salinities are a function of ground surface elevation and resulting inundation frequency only(5)Despite the young age of the Ostplate, a quasi-steady-state FWL has already been established. This is mainly a result of continuous mixing with intruding seawater due to the seasonal storm tides, which inhibits further growth of the FWL under the present geomorphological conditions at the narrow island(6)Real steady-state conditions will most likely never be reached due to the ever varying changes of geomorphology, sea level, and fresh recharge in such a highly dynamic coastal environment(7)A major uncertainty in the presented modeling approach are implemented recharge rates, especially the seawater infiltration rates. Future studies addressing the process of infiltration during inundation are needed

Data Availability

The observed and simulated data used to support the findings of this study are included in the article/Supplementary Material and are available from the corresponding author upon request, respectively.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

The authors wish to thank the Nationalparkverwaltung Niedersächsisches Wattenmeer (NLPV), especially G. Millat and C. Schulz for giving us the possibility to conduct studies on the Ostplate. Many thanks to the Niedersächsischer Landesbetrieb für Wasserwirtschaft, Küsten- und Naturschutz (NLWKN), mainly H. Dirks and R. Tants, for providing the laser-scan elevation models and the aerial image. We thank the Wasser- und Schifffahrtsverwaltung des Bundes (WSV) and the Bundesamt für Seeschifffahrt und Hydrographie (BSH) for providing data of the tide gauges. Many thanks are due to J. Alder, L. Diehl, T. Fresenborg, F. Grünenbaum, T. Haene, A. Harms, N. Kramer, C. H. Lünsdorf, H. Madaj, M. B. Maeso, S. Maiwald, J. Otten, F. Schindler, and S. Zahra for their support during field work. We would also like to thank U. Uebel and H. Freund for support during the construction of the observation wells. Many thanks to S. Fock, C. Heithecker, and C. Winkelmann of the Nationalparkhaus Wittbülten at Spiekeroog as well as to the team of the NLWKN located at Spiekeroog for their logistical support. Furthermore, we thank H. Nicolai and G. Behrens for the boat transfer. We also thank the German Research Foundation for project funding (DFG project number MA 3274/6-1). This publication is dedicated to our very supportive colleague Gerald Millat, who sadly passed away in 2017.

Supplementary Materials

The Supplementary Material consists of information regarding instrumentation and time of measurements (S.1), groundwater flow direction (S.2), grain size analysis (S.3), MHW-levels (S.4), temporal variability of the geomorphology at the Ostplate, recharge/infiltration rates, and recharge salinities implemented in the VG-model (S.5), present groundwater salinity distribution and extent of the FWL – Field data (S.6), videos of the temporal evolution of the groundwater salinity distribution and FWL for the CG-, VG-, and CGWST-model (S.7), comparison of the calculated and observed hydraulic heads (S.8), calculated groundwater salinities of the last five years of simulation (S.9), and the effect of choice of boundary conditions and parameter sensitivities (S.10). (Supplementary Materials)