Geofluids
Volume 2019 (2019), Article ID 1494326, 16 pages
https://doi.org/10.1155/2019/1494326

Research Article

Behaviour of Tritium and Tritiogenic Helium in Freshwater Lens Groundwater Systems: Insights from Langeoog Island, Germany

Vincent E. A. Post,¹ Georg J. Houben,¹ Leonard Stoeckl,¹ and Jürgen Sültenfuß²

¹Federal Institute for Geosciences and Natural Resources (BGR), Stilleweg 2, D-30655 Hannover, Germany
²University of Bremen, Institute of Environmental Physics, Bremen, Germany

Correspondence should be addressed to Vincent E. A. Post; vincent.post@bgr.de

Received 15 March 2019; Revised 13 July 2019; Accepted 21 July 2019; Published 7 October 2019

Academic Editor: Henrik Drake

Copyright © 2019 Vincent E. A. Post et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Tritium (³H) and its daughter product ³He have been widely used as tracers in hydrological studies, but quantitative analyses of their behaviour in freshwater lenses and the transition zone in coastal aquifers are presently lacking. In this paper, the fate of ³H and ³He in the freshwater lens and the transition zone as well as the saltwater wedge is studied using numerical variable-density flow and transport models of different degrees of complexity. The models are based on the conditions on the German island of Langeoog, which is uniquely suited for this purpose because of the high ³H concentration of the North Sea. It is found that most bomb-related tritiogenic ³He still resides in the freshwater lens, making it a useful tracer for young (<60 years) groundwater. Differences in dispersive transport between ³H and ³He can cause an apparent age bias on the order of 10 years. Under favourable conditions, ³H from seawater can penetrate deep into the offshore part of the aquifer and has potential to be used as a tracer to study saltwater circulation patterns. Our modelling suggests that the field-observed ³H in the transition zone does not originate from seawater but from freshwater affected by the bomb peak. Yet in models with a low () dispersivity, no ³H was sequestered into the transition zone and the transition zone width was underestimated. Better results were obtained with , a value that is higher than in comparable modelling studies, which suggests that further work is needed to better understand the controls (tides, lithological heterogeneity, or transience of recharge and pumping) on transition zone mixing processes.

1. Introduction

Tritium (³H) has been used extensively as a tracer in hydrogeological studies [1]. It is subject to beta-decay, and due to its half-life of 12.32 years, it is suitable to identify young groundwater. Tritium is naturally produced in the upper layers of the atmosphere but was also released in massive amounts by above-ground thermonuclear bomb testing in the 1950s and early 1960s [2]. Its atmospheric activity peaked around 1963-1964 and has been steadily falling since then.

When the tritium daughter product helium (³He) is simultaneously analysed and after the so-called tritiogenic contribution (³Hetrit) is separated from other ³He contributions (exchange with the atmosphere, excess air, and underground production), the apparent age (T) of a groundwater sample can be calculated based on the ratio of the ³Hetrit over the ³H concentration and the half-life constant:

in which is the radioactive decay rate of ³H. Strong variations in recharge rates over time and mixing of groundwater of different ages are a source of uncertainty for the calculated value of [3]. Lumped parameter models provide another means to determine the age [4, 5] and are particularly useful for dealing with mixing effects. They are most effectively applied when tracer time series are available.

While tritium has been used extensively in terrestrial, fresh groundwater settings, it has seen less application in brackish and saline coastal environments. It has been mostly applied in freshwater lenses [6–8] to calculate recharge rates with the Vogel [9] model:

where is the porosity, (L) is the lens thickness, and (L) is the depth. The growing interest in the use of brackish and saline groundwater as feed water for desalination plants [10], however, necessitates a reevaluation of the usefulness of isotopic tracers for studying saltwater dynamics.

Tritium has been detected in brackish and saltwater in coastal aquifers in Tahiti [11], Israel [12], Morocco [13], the Netherlands [14], China [15], and Australia [16]. Sivan et al. [12] inferred the timing of seawater intrusion by comparing extrapolated tritium decay curves for saline groundwater samples to the history of ³H concentrations in the Mediterranean Sea. Conclusive interpretations were not possible for the transition zone, as tritium may have derived from overlying freshwater and/or intruded seawater. That is, old, tritium-free seawater could have mixed with freshwater impacted by the bomb peak. In a similar manner, Han et al. [15] attributed the presence of ³H in brackish groundwater to the mixing of young groundwater with old, ³H-free saltwater from greater depth that is being drawn upwards by intensive groundwater abstraction.

While the tritium concentrations of most seas and oceans have returned to almost their prebomb era, natural background [17], a unique situation exists in northwestern Europe, where the discharge of water from nuclear processing plants sustains very high ³H levels in the North Sea [18]. As a result, Stuyfzand et al. [14] were able to attribute ³H values of up to ~30 TU (Tritium Units) in saline groundwater to the intrusion of ³H-rich seawater. The main motivation for the present work was the discovery of ³H in the transition zone between fresh and saline groundwater on the German island of Langeoog (Figure 1) during the study by Houben et al. [7]. They did not report these findings because their focus was on the freshwater lens, but tritium concentrations of up to 5.6 TU values were encountered in the transition zone that commences at 33 m below mean sea level (msl). It is conceivable that these values are due to the mixing of bomb peak-affected freshwater with intruded seawater, as put forward by Sivan et al. [12]. On the other hand, it may be that they originated from ³H-enriched seawater, similar to the findings by Stuyfzand et al. [14].

(a)

(b)

(c)

Figure 1: (a) Location of Langeoog Island, (b) map of the water supply area in the Heerenhus dune area showing the positions of the abstraction and observation wells, and (c) apparent bulk resistivity of the subsurface at 10 m below msl based on airborne electromagnetic measurements [22]. Blue lines in (b) represent the groundwater level elevation above msl during the summer of 1996, and the red line shows the projection of the model cross section. The shaded areas enclose the locations where Costabel et al. [23] found either no indication for the presence of the aquitard or found it to be <4 m thick. See their Figure 11 for the elevation of the top of the aquitard.

The relationships between tritium, helium, and groundwater age in coastal aquifers have not been explored using quantitative models except for the highly simplified case of the Henry Problem [19]. This lack of quantitative analysis means that interpretations of the age structure of freshwater lenses and intruded seawater bodies often remain ambiguous. The North Sea’s high tritium background concentrations provide an opportunity to study fresh and salt groundwater mixing processes along its coastlines. The objective of the present study is to extend the understanding of tritium dynamics in freshwater lenses in coastal aquifers. The focus is not only on the behaviour of ³H within the freshwater but also on the fate of ³H in the saltwater wedge and the offshore part of the aquifer. To this end, the spatiotemporal patterns arising from transient inputs and groundwater flow were investigated using numerical models, which are based on the hydrogeological conditions representative for Langeoog Island.

First, a set of models of increasing complexity was used to determine the influence for each of the various controlling factors separately. Attention is paid to the implications for submarine groundwater studies. The results are expected to be applicable to other study areas, especially along the North Sea shores, and also to other coastal localities where tritium concentrations persist above the natural background (e.g., [20]). Second, following a discussion of the model results, we address the issue of the origin of ³H in the transition zone (i.e., whether it originates from the freshwater or the seawater). In what follows, the notation ³He is used to indicate the tritiogenic ³He contribution, omitting the subscript trit.

2. Study Area

The island of Langeoog is part of the chain of barrier islands that separate the shallow intertidal Wadden Sea from the North Sea (Figure 1). It first formed around 2,800 to 2,200 years before present as a sand bar, onto which dunes developed around the year 100 AD [21]. These remained largely unvegetated due to constant wind erosion until the 13th century. Thereafter, the formation of soils allowed the first human settlements. Over the following centuries, dune stabilization by vegetation and earthworks led to the formation of two parallel, west-east trending dune belts, which continue to dominate the island morphology.

Below the superficial dune and beach sands, the subsurface is made up of Holocene deposits of the Wadden Sea, underlain by mostly Pleistocene (glacio-) fluvial sediments [24]. Intercalated low-permeability sediments form an aquitard, which, because of glaciotectonic activity, is known to be spatially discontinuous (Figure 1(b)) and highly variable in thickness. In the northern part of the study area, the low-permeability unit consists of clay. No information is available about its offshore continuation. In the southern part, the unit is more heterogeneous and consists of fine sand, silt, and clay. It generally occurs below 20 m below mean sea level and, where well developed, it can reach a thickness of 10 m or more [23].

Three individual freshwater lenses exist on the island. Their separation is the result of occasional catastrophic storm flood inundations, storm floods, e.g., the Christmas flood of 1717, which caused large breaches of the dune belt in two locations [25]. Only the western lens is currently used, the outline of which can be inferred from the measured bulk apparent resistivity (Figure 1(c)) obtained from helicopter airborne electromagnetic (HEM) surveys [22]. Before the rise of tourism, the population obtained its water from rainwater and shallow dug wells. To satisfy the increasing demand, the first three pumping wells were drilled in 1909 close to the village. After 1938, they were augmented and later replaced by wells drilled east of the village in the Heerenhus dunes, where a well field is still active today (Figure 1). Most wells are located in the Pirola Valley, which is situated between the two dune belts stretching W-E along the island’s northern shoreline.

Extraction is done intermittently at pumping rates of 10 m³/hr per well and is distributed over 20 small wells, with screens installed at varying depths between 10 and 18 m below sea level. As an effect of increasing tourism, water consumption continuously rose until the 1980s. It has significantly decreased since the 1990s due to the implementation of water-saving measures. Peak demand occurred in 1983 with 452,000 m³/yr (Figure 2), while by 2011 the pumped volume was down to 333,000 m³/yr, despite a mostly stable visitor number to the island since 1990 [7].

Figure 2: Time series data used in the model: (a) total annual abstraction () from the freshwater lens (blue bars) as well as seasonally distributed abstraction per unit of surface area () versus time (red line); (b) calculated groundwater recharge (, bars) and ³H concentration of groundwater recharge (³Hrech, solid red line) and seawater (³Hsea, dashed red line) as a function of time (see text for source of data). The colour of the bars in (b) indicates relatively dry (yellow) or wet (green) phases based on the cumulative deviation from the mean recharge rate. The small green arrow pointing to the axis in (b) indicates the mean recharge rate () that was used in simulation A1.

3. Methods

3.1. Sampling and Analysis

Basic information, such as well locations, groundwater levels, and pumping rates, were obtained from the local water supply company, the Oldenburgisch-Ostfriesischer Wasserverband (OOWV). Details about sampling and analytical techniques used to determine the chemical and isotopic composition of the groundwater samples were presented by Houben et al. [7]. He isotopes and tritium were analysed at the Institute of Environmental Physics of the University of Bremen [26]. The low detection limit for tritium (0.02 TU; ratio of ) was achieved by using the ³He accumulation method. For all sampling sites, duplicate samples were taken and analysed.

Transition zone water samples were obtained from observation well 17 (Figure 1), which has a long screen (from 33.6 to 57.6 m below sea level) reaching into the brackish water. Profiles of specific conductance (, i.e., the electrical conductivity for a reference temperature of 25°C) versus depth were obtained on two occasions, using an Ott™ Hydrolab multiparameter probe. After the first measurement, on 30 May 2013, the well was purged. An almost identical profile was afterwards established, as evidenced by a second measurement on 15 November 2013 (shown in Figure 3(a), alongside with model results that will be discussed later on). During the second sampling, three individual water samples were obtained by low-flow sampling using a Grundfos MP1 pump from three different depths, similar to Lee et al. [27]. It is known that coastal zone observation wells with long screens may indicate a higher position of the transition zone than in the aquifer [28, 29]. The resemblance between the specific conductance values of the three samples and the probe-derived values, as well as the HEM data [23], provides confidence though that the borehole data are representative for the groundwater.

Figure 3: (a) Chloride concentration, (b, c) ³H concentration, and (d) ³He concentration as a function of depth for observation well 17. The small red square marker symbols in (a) indicate the chloride concentration calculated from the downhole log. Round symbol markers in (a), (c), and (d) represent samples collected in 2013 (previously published by Houben et al. [7]); square symbol markers are for the samples from well 17. Crosshair symbols in (b) are data for January 2002 from Marggraf [30]. The displayed samples are from different wells within 300 m of the model cross section. The solid lines represent the model outcomes (discussed later on) of simulations A7-A9; the dashed lines are for simulation A7 with longitudinal dispersivity . See Figure 1(b) for the location of observation well 17 and Figure 10 for its projected location on the modelled cross section.

3.2. Numerical Modelling

A two-dimensional section across the freshwater lens was considered (Figure 4). SEAWAT [31] was used to solve the following mass balance equations for fluid and solute, respectively,

where is the fluid density (M L⁻³), is the Darcy flux (L T⁻¹), is the density of water entering (or leaving) through a source (or sink) [M L⁻³], is the source’s volumetric flow rate per unit volume of aquifer (T⁻¹), is the concentration (M L⁻³), is the hydrodynamic dispersion tensor (L² T⁻¹), and is the solute concentration associated with water sources and sinks (M L⁻³). The hydrodynamic dispersion tensor is a function of the longitudinal () and transverse () dispersivity (L), as well as the effective molecular diffusion coefficient () (L² T⁻¹).

Figure 4: Cross section showing the model dimensions, boundary conditions, and hydraulic parameters for the layers used in the layered simulations. Parameter values in white font were constant in all model simulations. Along the left, bottom, and right boundaries, no-flow and nondispersive flux conditions were applied. The abstraction () and recharge () time series used are shown in Figures 2(a) and 2(b), respectively. The input time series for ³H in recharge (³Hrech) and seawater (³Hsea) are shown in Figure 2(b). The fading white colour of the clay layer indicates the uncertainty about its offshore extent.

The term (M L⁻³ T⁻¹) in equation (4) accounts for the effects of radioactive decay. Three dissolved species, Cl, ³H, and ³He, were simulated. For Cl, , while for ³H, the decay rate is , where the subscript H is used to indicate that C represents the ³H concentration. For ³He, , i.e., the growth rate of tritiogenic ³He is equal in magnitude but opposite in sign to the decay rate of ³H. A modified version of the SEAWAT code was used to simulate this type of parent-daughter chain reaction based on the same approach as Bedekar et al. [32].

Nine different scenarios (A1-A9) were considered that each differ from each other by a single model variable (Table 1). The simplest simulation (simulation A1) had a constant recharge rate, a homogeneous and isotropic hydraulic conductivity, and no groundwater abstraction. The complexity of the model simulations A2 to A9 was increased in a stepwise manner as shown in Table 1. The most complex simulations (A7-A9) resemble the conditions on Langeoog Island the closest. In A8, the clay layer was made continuous in the offshore domain, but in A9, it was truncated 200 m offshore as there is considerable uncertainty about the spatial extent and the hydrodynamic significance of this layer [23]. The choice for 200 m was therefore also simply arbitrary. The model parameters and the boundary conditions are shown in Figure 1. A set of trial simulations was conducted to determine the model thickness at which the freshwater lens and the transition zone were not affected by boundary effects. The thickness of the aquifer sublayers was determined based on available borehole data [21, 23]. The hydraulic conductivities were taken from a calibrated three-dimensional finite element model for the entire island by Stoeckl [33], who used hydraulic conductivity estimates based on particle size by Marggraf and Naumann [34] and unpublished pumping test data.

Table 1: Overview of the model simulations.

A1
A2 (A3₁)
A3
A4
A5
A6
A7
A8
A9

Transient recharge

X
X
X
X
X
X
X
X

Anisotropic

X
X
X
X
X
X
X

Layered

X
X
X
X
X
X

Clay layer

X

X

Discontinuous clay layer

X

X

Abstraction

X
X
X

The effect of the vertical anisotropy on the ³H concentration distribution in the offshore domain due to the recirculation of seawater was investigated in an additional set of simulations labelled A3. The horizontal to vertical hydraulic conductivity ratio was varied between values . This number was used as a subscript to indicate the simulation, e.g., A3₅ is for simulation A3 with . Since simulation A2 differs only from A3 in that it has , it has been equivalently labelled as A3₁. Some simulations with a constant ³H input or a higher dispersivity were also run but not specifically labelled.

The transient simulation spanned the 61-year period from 1 January 1953, the earliest month for which the rainfall ³H concentration could be estimated, to 31 December 2013. Groundwater recharge was calculated based on daily rainfall and evaporation measurements as described in Post and Houben [35], and monthly averages were used in the numerical model (Figure 2(b)). The ³H concentration of the recharge was estimated from the measured concentrations of nearby rainfall stations of the global network of isotopes in precipitation (GNIP) [36]. The partial records of the German stations in Stuttgart and Cuxhaven as well as the Dutch station in Groningen were collated to form a continuous input time series (Figure 2(b)). Being the nearest to Langeoog (about 82 km), the Cuxhaven data (available from January 1978) were used as the principal time series. The time series was further completed with the Groningen data (85 km away, available from January 1970) and those from Stuttgart (available from February 1962). For the period prior to 1962, the data from Ottawa, Canada, were used. Due to the usually small depth to groundwater, it was assumed that the travel time across the unsaturated zone was short, so that the rainfall tritium concentration for a particular month could be assigned to the recharge of that month.

The ³H concentration of seawater (Figure 2(b)) was estimated from different sources. Values prior to 1977 were assumed to be equal to those measured in the Atlantic Ocean [37, 38], while after 1995 a constant value of 30 TU was assumed, which is due to the disposal of tritium-enriched water into the English Channel from nuclear power and processing plants at the French coast [18]. This water flows from the English Channel along the southern coast of the North Sea to the Norwegian coast in a counter-clockwise pattern. A linear increase was assumed between 1977 and 1995, which was found to be in agreement with ³H concentrations measured in the North Sea [39, 40] that can be attributed to increased emissions from nuclear power plants in France, Belgium, and Germany into rivers feeding the North Sea.

The initial conditions for the model (i.e., the heads and concentrations on 1 January 1953) were determined by simulating the 61-year period twice (i.e., 122 years spanning the period 1 January 1831–31 December 1952). With the exception of simulation A1, the initial conditions therefore do not represent a steady-state situation, the reason being that the temporal variability of the recharge could not be neglected as it is responsible for a significant widening of the transition zone. Using a constant recharge to generate a steady-state starting point would thus result in an unrealistically narrow transition zone, as discussed below. During the 122-year period, all model parameters remained the same except that no abstraction was considered and that the ³H concentration of the meteoric recharge was assigned a constant value of 5 TU, the estimated natural background (Craig and Lal, 1966; Roether, 1967). It was further assumed that seawater had prior to 1953 [38].

To include abstraction in the two-dimensional cross-sectional simulations, the reported abstraction volumes were divided by the estimated ground surface area of the well field zone of influence of . Groundwater withdrawal was distributed over the model cells between 200 and 400 m from the left boundary and 9.5 to 14.5 m below mean sea level (a total of 75 cells). This approach assumes that the effect of the withdrawal extends uniformly in the direction perpendicular to the cross section, which is a justifiable approximation as the water table elevation shows little variation in that direction (Figure 1(b)). Because pumping is highly seasonal, but the exact monthly amounts were not available for most of the years, an estimate of the intra-annual distribution of pumping rates was made based on monthly figures for the years 1995, 2001, and 2011. Based on these data, pumping was found to be lowest in January and February (4% of the annual total) and highest during the tourist season in July and August (14% of the annual total). Given the constancy of both climate and tourism seasons, it could be safely assumed that this distribution was the same for all the years [7], so monthly abstraction rates were inferred from the annual totals for the entire simulation period in this way (Figure 2(b)).

Post and Houben [35] successfully modelled downward descending salt fingers on the neighbouring island of Baltrum using , but since their model had a much smaller scale, it was considered to be too small for the present study. The longitudinal dispersivity was therefore set to (with ). This value is intermediate between the values used by Stuyfzand [41] and Vandenbohede and Lebbe [42], who adopted values of and , respectively, for the freshwater lenses in the Dutch and Belgian coastal dunes, and the value used by Pauw et al. [43] who used for the larger barrier island of Texel in the Dutch part of the Wadden Sea. The effect of this parameter will be discussed in more detail below. Effective molecular diffusion coefficients of , , and were used for Cl, ³H, and ³He, respectively [44].

4. Results and Discussion

4.1. Transient Effects

The results of simulations A1 and A2 are shown in Figure 4. It can clearly be seen that the temporal variability of the recharge results in a widening of the transition zone between fresh and saline groundwater in simulation A2 (compare Figures 5(a) and 5(b)). At the left model boundary, the distance between the chloride contour lines representing 2.5 and 97.5% of the seawater chloride concentration is less than 8 metres for simulation A1, whereas it is ~25 metres for simulation A2. The variability of the recharge causes a net groundwater movement that promotes mixing between the fresh- and saltwater. This effect is similar to the widening of the transition zone caused by tidal fluctuations [45].

This result highlights the importance of explicitly modelling rainfall recharge as a transient process for understanding the development of the transition zone in coastal aquifers. The importance of the temporal rainfall recharge variability has also been recognized in other modelling studies of island aquifers. For example, Griggs and Peterson [46] modelled the freshwater lens of the Laura area on the Majuro Atoll, Marshall Islands. They found that the high temporal variability of the recharge rate resulted in an absence of steady-state conditions and concluded that the sustainable yield could not be modelled if yearly instead of monthly averages were used. Both Ghassemi et al. [47] and Post et al. [48] were only able to reproduce the temporal variability of the salinity of abstracted groundwater from an atoll freshwater lens if they adopted monthly averaged recharge rates.

The effect of the ³H bomb peak can clearly be recognized as a quasihorizontal, elongated zone of high ³H concentrations, which starts at a depth of 20 to 40 m below sea level near the left model boundary and gradually rises closer to the surface with decreasing distance to the shoreline (Figure 5). This pattern is due to the freshwater flow pattern within the lens. The peak concentrations in simulation A2 are lower than those in A1, which is partially because of the dispersive spreading effect caused by the temporally variable recharge rates. Another important reason though is that around the time when the atmospheric ³H concentrations were at their highest, groundwater recharge rates were relatively low (Figure 2(b)). Therefore, the total mass of ³H that entered the freshwater lens was considerably lower and resulted in a lower total ³H mass value for simulation A2 than for simulation A1 (Figure 6). For example, at the peak in September 1966, total ³H contained within the freshwater lens amounted to for simulation A2, compared to for simulation A1.

(a)

(b)

(c)

Figure 5: Cross sections showing contour plots of the simulated ³H concentrations (shaded colours) and Cl concentrations (white lines) at the end of the simulation period (31 December 2013). The five white lines represent Cl concentrations equal to 2.5%, 16%, 50%, 84%, and 97.5% of the seawater chloride concentration (17 g/L). Results are shown for simulations (a) A1, (b) A2 (=A3₁), and (c) A3₂₀. The dashed white line in (a) and (b) indicates the position of the vertical profiles shown in Figure 7.

Figure 6: Temporal evolution of the total ³H in Bq (solid lines) and ³He (dashed lines) in the freshwater lens (groundwater with less than 2.5% seawater) for simulations A1 (red) and A2 (blue).

The graph in Figure 6 further shows the total ³He produced from tritium in the freshwater lens. The maximum was reached during the mid-1980s but, as ³He production in the freshwater lens fell because of the declining ³H concentrations, the ³He lost by discharging groundwater started to exceed production and total ³He has been decreasing since then. Nevertheless, most (about 80%) of this ³He was still present within the freshwater lens by the end of the simulation. Groundwater dating back to the 1960s thus remains identifiable by high ³He concentrations it acquired due to the ³H peak in recharge at that time. This provides opportunities to use ³He in lieu of ³H to characterise groundwater flow patterns in freshwater lenses where the groundwater residence times are on the order of a few decades.

Solomon and Sudicky [49] investigated the effect of ³H in rainfall being variable in time on the apparent age . Tritium enters the saturated zone with the recharge, but the decay product ³He forms in situ in the aquifer. When the input signal of ³H is steady, steady-state ³He and ³H concentration profiles form, and the effect of dispersive mixing of groundwater of different ages on the calculated value of is small [49]. When the input of ³H is variable in time, such as in the case of the ³H bomb peak, ³He is being produced at different rates in different parts of the aquifer. The resulting concentration gradients drive dispersive and diffusive transport away from the zones of the highest ³He concentrations. As the concentration gradients of both ³H and ³He vary with time, this results in temporally variable dispersive and diffusive transport rates of both tracers. Solomon and Sudicky [49] found that this can result in an error of of up to 5 years in the groundwater above the tritium peak. Running simulations A1 and A2 with a constant atmospheric background ³H input (simulations A1c and A2c) and comparing the differences show that groundwater above the peak is seemingly older by 10.3 and 8.7 years (simulations A1 and A2, respectively) (Figure 7). The model parameters used in this study are comparable to those used by Solomon and Sudicky [49], and the higher age mismatch found here is therefore mainly because in our case the comparison is made about 50 years after the bomb peak, versus about 20 years in their study. Solomon and Sudicky [49] also found younger apparent ages below the tritium peak, and this is also the case here (Figure 7). Because groundwater ³H concentrations in this part of the lens have fallen to below detection limits, the difference is no longer of any practical consequence. This bomb peak-impacted groundwater remains identifiable though based on its elevated ³He concentration.

Figure 7: (a) Apparent age (solid lines) and chloride concentration (dashed lines) for simulations A1 and A2 and (b) the ³H (solid lines) and ³He (dashed lines) concentration versus elevation at the end of the simulation period (31 December 2013). The location of this profile is indicated in Figure 5. Results are shown both for a constant ³H input of 9 TU (as indicated by the subscript c) and for the ³H input according to the curve in Figure 2(b).

Comparing simulations A1c and A2c shows that the transience of the recharge similarly results in a small error of the calculated ³He/³H ages even when the ³H input signal is constant in time. In the case of simulation A2c, the maximum difference of the ³He/³H age compared to simulation A1c is 1.8 years. The difference is due to the greater dispersion in simulation A2c than in simulation A1c, but the magnitude of the effect is only of secondary order compared to that of the bomb peak ³H input.

The apparent age shift due to the nonconstant ³H input is also apparent when the sum of ³H and ³He is plotted as a function of (equation (1)), which is often done to identify mixing effects in field data (e.g., [8]). If ³H and its daughter ³He were to remain within the same parcel of groundwater, their sum would reflect the ³H concentration at the time of recharge. In that case, a groundwater sample would have a matching data point somewhere on the atmospheric ³H concentration input signal in the graph of Figure 8. In reality of course, dispersive mixing causes attenuation of the peak input concentration. Moreover, because of the aforementioned dispersive separation of ³H and ³He and its effect on , the peak of the summed modelled ³H and ³He concentrations of the groundwater does not coincide with the peak of the atmospheric ³H concentrations. Both simulations A1 and A2 fail to match the measured data of samples recharged before 1974, which could be because the models underestimate mixing effects, which will be discussed further below.

Figure 8: Comparison of the ³H concentration in precipitation with time (same data as shown in Figure 2(b)) to the sum of the simulated and measured ³H and ³He concentrations in groundwater versus for simulations A1 (red dots), A2 (blue dots), and A7 with (purple dots). Each dot represents the simulated value for a model cell. Only values of cells or samples with less than 2.5% seawater are shown.

4.2. Tritium Behaviour in the Coastal Discharge Zone

The recirculation of seawater results in an asymmetric funnel-shaped zone of elevated ³H concentrations in the saline groundwater just seaward of the coastline (Figure 5). In simulation A1, this zone even reaches below the transition zone at the bottom of the freshwater lens. In simulation A2, the ³H-enriched seawater penetrates to an even greater depth, albeit that it does not reach below the transition zone anymore, because the latter expands more strongly. Since the circulatory flow of seawater [50] is driven by dispersive mixing [51, 52], the reconfiguration of the zone of high ³H concentrations is expected and can be attributed to the enhanced mixing in simulation A2 compared to simulation A1. These results also emphasise the difficulty of using lateral changes of ³H (or ) to infer seawater intrusion rates, because of the strong vertical age stratification within the saltwater wedge [53]. This complicates the selection of a saltwater end member to correct for mixing effects in the interpretation of age tracers in coastal aquifers [54, 55] and calculating the rate of seawater intrusion [12].

An important factor that determines the circulatory flow pattern of seawater is the vertical hydraulic conductivity [52]. The effect of this parameter was explored in a set of models based on simulations A2 and A3 by letting the anisotropy ratio be 1, 2, 5, 10, and 20 by varying . It is noted that the simulation results show that the freshwater lens is affected only to a small degree (compare Figures 5(b) and 5(c)), because there is little change in the groundwater pressure distribution within the lens between the simulations. The reason for this is that the latter is primarily controlled by and the recharge , which both remained unchanged between the simulations, and consequently, the water table elevation changed by less than a centimetre at the inland boundary. Near the coast, the differences were larger, and this is reflected by a different shape of the transition zone in the interval .

In contrast to the relative constancy of the freshwater lens geometry, Figure 5(c) clearly shows that the decrease of significantly reduces the maximum penetration depth of the funnel-shaped body of high ³H saltwater (measured by , the elevation of the deepest point where TU) below the seafloor. This indicates an adjustment of the flow field in the subsea part of the aquifer. To analyse this in more detail, the flux across the top model boundary in the specified-head cells (representing the seafloor) was calculated. As the magnitude and direction of the flux depend strongly on the recharge rate [48, 56], the fluxes at the end of each monthly stress period were averaged over the duration of entire simulation. Discharge of land-derived freshwater mixed with recirculated seawater occurs across the cells nearest to the coast. The width of the discharge zone widened from one cell (10 m) to 7 cells (70 m) as changed from 1 to 20. The effect on the inflow of seawater across the seafloor () is shown in Figure 9. The shape of the relationship between and is the same as the relationship between and defined by Smith [52] (his Figure 2), where . In the setup used by Smith [52], and represented the lateral freshwater inflow across the inland model boundary per unit shoreline length and per unit cross-sectional area, respectively. In our simulations, the freshwater derives from recharge across the top model boundary. Nevertheless, the resemblance of the results is expected, as only changed between the simulations and all other parameters that determine remained constant.

Figure 9: Dependency of (a) seawater inflow and (b) (the elevation of the deepest point where ³H> 0.02 TU) on vertical hydraulic conductivity for simulation A3 with different anisotropy ratios. The corresponding anisotropy ratios are displayed along the top horizontal axis of the graph. was calculated by summing the inflow across the specified head cells along the seafloor. The vertical error bars represent one standard deviation of the temporal variability of and , which is caused by the transient recharge (Figure 2(c)).

The practical implication of this finding is that for studies of submarine groundwater discharge it may be possible to use ³H to constrain the recirculated seawater component. For conditions similar to those simulated here (i.e., a sufficiently homogeneous aquifer) and assuming the seawater is well mixed and the vertical anisotropy is constrained, ³H measured in water samples from a multilevel observation well at the shoreline might be used to infer the magnitude of . A relatively high ³H content of seawater, as in the case of the North Sea, is clearly advantageous for that purpose.

4.3. Heterogeneity Effects

In simulation A4, the horizontal hydraulic conductivity of the upper 40 m was lowered (Figure 3), in accordance with the available geological information for Langeoog. Predictably, this resulted in a significant thickening of the lens (Figure 10(a)) which is caused by an increased water table elevation (i.e., higher groundwater pressures) as the freshwater in the lens experiences a greater resistance as it flows towards the sea. The presence of a continuous or discontinuous clay layer (Figures 10(b) and 10(c), respectively) adds to the resistance against freshwater flow to the sea and thereby caused the freshwater lens to thicken by a few metres in simulations A5 and A6 relative to simulation A4. The higher conductivity layer beneath 40 m depth caused a minor, barely perceptible dissection of the main ³H plume, as flow velocities in this layer are higher than in the overlying less permeable layer. In both simulations A5 and A6, the freshwater lens and the transition zone extended into the offshore part of the aquifer, which is attributed to a diversion of the freshwater flowing towards the coast. The clay layer prevents it from flowing upward, and instead, a part of the freshwater is forced to flow laterally underneath the clay. In simulation A5, the freshwater tongue beneath the clay layer had not reached a stable position by the end of the simulation and kept slowly expanding in a seaward direction.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10: Contour plots of the simulated ³H concentration activities (shaded colours) and Cl concentrations (white lines) at the end of the simulation period (31 December 2013). The five white lines represent Cl concentrations equal to 2.5%, 16%, 50%, 84%, and 97.5% of the seawater chloride concentration (17 g/L). Results are shown for simulations (a) A4, (b) A5, (c) A6, (e) A7, (f) A8, and (g) A9. The location of the clay layer is indicated by the hatched white rectangle. The open red rectangle in (d)–(f) encloses the cells with groundwater abstraction, and the dashed white line indicates the projected location of observation well 17.

A complex distribution of chloride and ³H developed above the clay layer. The reason for this is that the clay has a dispersing effect on the groundwater discharge near the coast. Rather than being focussed in a relatively narrow zone, the discharge extends over circa 200 m in simulations A5 and A6 (as well as A8 and A9). Hence, the upward flow rates decrease, and the propensity for density-driven downward flow becomes greater. Consequently, a complex pattern of upwelling and downwelling occurred, and water fluxes across the aquifer-seawater interface were spatially and temporally variable. Similar effects have been noted in numerical simulations of submarine groundwater discharge by Kooi and Groen [57]. While the complex dynamics of convective fingering may not be accurately captured by numerical models at this scale, this result indicates that some of the variability of SGD may be due to the presence of large-scale geological features and that small-scale heterogeneity is not necessarily a precondition for spatial variations of the SGD flux. Similar behaviour has also been observed in physical sand tank experiments of heterogeneous coastal aquifers [58]. When these conditions occur, the presence of ³H below the seafloor in the discharge zone is associated with downwelling salt fingers. The funnel shape observed in simulations A1–A4 only developed seaward of the discharge zone, although in simulation A5 (as well as A8) its development was very much suppressed due to the presence of the continuous clay layer (Figure 9(b)).

With respect to the role of the clay layer, it is also interesting to note that simulations A5 and A6 had almost the same chloride and ³H concentration patterns below the dune area (Figure 10). This means that the uncertainty about the offshore geology did not impart great uncertainty on the onshore solute transport simulation results. This is attributed to the fact that almost stagnant flow conditions exist in the aquifer below the seafloor, and therefore, the presence of a low-permeability unit did not alter the flow field significantly. Given that there is a lack of offshore geological data, this is a favourable outcome, but to what extent this holds for coastal aquifer studies more generally is uncertain.

Simulations A7–A9 included the effect of groundwater abstraction. Comparing Figures 10(a)–10(c) to 10(d)–10(f) shows that the effect was that the freshwater lens thinned by about 6 m, as measured by the vertical position of the 2.5% seawater salinity contour at the left model boundary. The distribution of ³H was affected in two ways: first, the upward movement below the zone of pumping drew in water with bomb peak tritium from the deeper parts of the lens. Second, ³H activities were generally lower in the simulations with abstraction. This is attributed partially to the removal of water containing ³H through the abstraction wells, but it cannot be excluded that enhanced dispersive spreading due to the more irregular plume shape may have played a role. A widening of the fresh-saltwater transition zone below the abstraction wells can be recognized, which is consistent with observations in laboratory sand tank experiments of pumping-induced up-coning [59].

4.4. Comparison to Measured Data

The simulated and measured variations of the chloride concentration as well as the ³H and ³He concentrations with depth are shown in Figures 3(a)–3(d), alongside with the data published by Marggraf [30] and Houben et al. [7]. Arguably, a direct comparison between the simulation results and field data is difficult because of the limitations of the two-dimensional model and the required projection of the samples onto the modelled transect, but nonetheless, the simulations’ resemblance of the overall observed trends with depth is encouraging. Visually, simulation A7 appears to better match the ³H data in the upper 20 m of the profile as measured in 2013 than simulations A8 and A9 (Figure 3(c)). Simulation A7 also appears to best capture the ³H peak observed in January 2002 (Figure 3(b), note the difference in scale with the graph of Figure 3(c)) as well as the measured ³He activities (Figure 3(d)).

This outcome is somewhat surprising as simulation A7 did not include the clay layer, which is known to exist in the area albeit that there is uncertainty about its hydrogeological significance. Possibly, this result indicates that simulations A8 and A9 overestimated the effect of the clay layer on solute transport in the lens. This might be because its permeability is higher than the value used in the model (Figure 4) or that the clay layer is too discontinuous to act as a real aquitard. It is known that its continuity has been affected by glaciotectonic effects [23]. Further simulations, based on a three-dimensional representation of the system, are required to better understand the role of the clay layer and the effect of abstraction.

The depth of the transition zone at the projected location of well 17 was overestimated in simulations A8 and A9. The midpoint of the transition zone was well matched by simulation A7 (Figure 10(e)), but the simulated width was too small relative to its observed width. It is likely that the field data are not an artefact caused by the long well screen or sampling procedures, because geophysical data (HEM and time-domain electromagnetic measurements) for Langeoog Island also suggest the presence of a thick (20 to 30 m) transition zone beneath the freshwater lens [23]. This implies that dispersive mixing between the freshwater and saltwater was underestimated by the model. Simulation A7 was therefore run with a longitudinal dispersivity (i.e., ten times the original value, and increased by the same factor). The results in Figure 3(a) show that the simulated transition zone width closely resembled the observed width (Figure 3(a)). While the peaks for ³H and ³He in the overlying freshwater lens became more attenuated, the model results still resembled the patterns displayed by the field data.

In none of the simulations did ³H originating from intruded seawater reach the location of the observation well, which most likely rules it out as a source of the observed ³H in the transition zone at appreciable inland distances from the shore. This contrasts the findings by Stuyfzand et al. [14], but in their case, a significant horizontal inland movement of the saltwater wedge occurred by excessive pumping, which is not seen in our models. The model with a longitudinal dispersivity did result in elevated ³H concentration in the transition zone, suggesting that the mechanism put forward by Sivan et al. [12], i.e., freshwater containing bomb-derived ³H mixes with ³H-free seawater, can be responsible for ³H in the transition zone in coastal aquifers.

But even with a longitudinal dispersivity , the Langeoog models underestimated the measured ³H concentrations in the transition zone. The underrepresentation of mixing processes is most likely because important confounding factors were not simulated. The consideration of tides is known to lead to a significant widening of the transition zone [45] but is computationally too expensive at the decadal timescale of our models. Multidomain solute mass exchange effects can also cause a widening of the transition zone [60] and may be important here as the tidal deposits are characterised by decimetre scale grain size variations. Also, even though pumping rates were varied using monthly time steps (Figure 2(a)), the temporal dynamics caused by pumping may be even greater; as in practice, individual wells are operated on a rotational basis whereby they are operated for periods lasting up to 24 hours and are then switched off. This is done to minimise up-coning, but the highly transient flow field that results from this way of operating the well field may induce greater mixing than what is being simulated in the model. Moreover, the combination of three-dimensional flow and subsurface heterogeneity may lead to more dispersion than was captured by the two-dimensional simulations. Dispersivity values greater than were not considered but may be required to explain the width of the transition zone.

5. Conclusions

The model results presented in this paper demonstrate how atmospherically derived tritium behaves in freshwater lens groundwater systems under different combinations of recharge variability, aquifer anisotropy, lithological heterogeneity, and groundwater abstraction. The results reinforce earlier findings by Solomon and Sudicky [49] that the highly transient input of bomb ³H in groundwater recharge results in a significant age bias for groundwater recharged around the time of the bomb peak in the 1960s, which is due to the difference in dispersive spreading of ³H and its daughter ³He in the aquifer.

The model results can be synthesised in terms of groundwater age classes as shown in Table 2, which is expected to be applicable to other freshwater lens systems along the shorelines of the southern North Sea. A finding that is expected to be generally applicable is that, since most of the tritiogenic ³He was found to still be present within the freshwater lens, it may be possible to use ³He as a tracer to identify freshwater recharged around the time of the bomb peak (a role traditionally fulfilled by ³H). Because ³He was also sequestered into the transition zone, it can also be used to study fresh-saltwater mixing relationships, which can contribute to the understanding of the formation of brackish groundwater in coastal areas.

Table 2: Time since recharge for freshwater and seawater depending on groundwater ³H and ³He concentration (all numbers in TU) for simulation A7 with .

Freshwater
≥60 yr (³He ≈ 5)
<10 yr
≥20 yr
<10 yr ()

Saltwater
≥40 yr (³He ≈ 0.2)
Mixture of old seawater and
Mixture of seawater and ≥20 yr freshwater (brackish)
<10 yr ()

Using a transient recharge input as opposed to a constant recharge had only a secondary effect on apparent ages but did result in a significant widening of the transition zone between fresh- and saltwater. While the width of the transition zone (as measured at well 17) could be approximated by increasing the dispersivity by an order of magnitude (simulation A7), the simulated ³H concentrations for the transition zone underestimated the measured values. This suggests that the model still underrepresented mixing between fresh- and saltwater, despite the adopted dispersivity value of being higher than what was used in modelling studies of comparable freshwater lenses [41–43]. The limitations of the present model complicate a direct comparison to the field data, and a fully three-dimensional model with appropriate consideration of local heterogeneity effects is required to better understand the formation of the wide transition zone. If time series data become available, lumped parameter models could further be used to investigate the role of dispersion [61].

Despite uncertainty about the magnitude of the dispersion coefficient and processes that control the width of the transition zone, our models suggest that the ³H observed in the transition zone on Langeoog more likely derived from the freshwater lens than from the North Sea. In none of the simulations did ³H from the North Sea reach far enough inland to explain the observed values.

In the offshore part of the aquifer, there can be a pronounced funnel-shaped zone of elevated ³H concentrations that stems from the circulation of seawater driven by dispersion in the transition zone. The vertical hydraulic conductivity and the presence of a clay layer were found to exert a strong control on the maximum depth to which measurable ³H can be found in the subsea portion of the aquifer, as well as on submarine groundwater discharge patterns. Therefore, groundwater ³H measurements at the coastline or below the seafloor could see application as an indicator of the strength of the seawater circulation and SGD in the southern North Sea area and other coastal sites where ³H or ³He variations can be resolved with sufficient analytical precision.

Data Availability

Previously reported field data were used to support this study and are available at [doi 10.1002/2014WR015584, doi: 10.1016/j.jappgeo.2016.11.007, and https://d-nb.info/975110349/34]. These prior studies (and datasets) are cited at relevant places within the text as references [[7, 23]; and [30]].

Additional Points

Key Points. Numerical models provide insight into the behaviour of ³H and ³He in island aquifers. Bomb-related tritiogenic ³He still provides a dating tool even though all ³H has decayed. Bomb-pulse ³H input causes an apparent age bias on the order of up to 10 years. Saltwater ³H distribution is sensitive to seawater recirculation and submarine groundwater discharge. Sequestration of bomb-pulse ³H in transition zone requires vigorous mixing.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to thank the Landesamt für Bergbau, Energie und Geologie (LBEG), of Lower Saxony for the technical support during the sampling of the wells on Langeoog and the Oldenburgisch-Ostfriesischer Wasserverband (OOWV) for providing data. Dr. Matthew Currell provided thoughtful reviewer comments that led to significant improvements of the manuscript.

Supplementary Materials

This file contains a table with the data previously published by Houben et al. [7], extended with the data for well 17 in the transition zone. (Supplementary Materials)

References

D. K. Solomon and P. G. Cook, ³H and ³He. Environmental Tracers in Subsurface Hydrology, P. G. Cook and A. L. Herczeg, Eds., Springer US, Boston, MA, USA, 2000.

F. Eyrolle, L. Ducros, S. Le Dizes et al., “An updated review on tritium in the environment,” Journal of Environmental Radioactivity, vol. 181, pp. 128–137, 2018.

B. C. Jurgens, J. K. Böhlke, L. J. Kauffman, K. Belitz, and B. K. Esser, “A partial exponential lumped parameter model to evaluate groundwater age distributions and nitrate trends in long-screened wells,” Journal of Hydrology, vol. 543, pp. 109–126, 2016.

B. C. Jurgens, J. K. Böhlke, and S. M. Eberts, TracerLPM (version 1): an Excel® workbook for interpreting groundwater age distributions from environmental tracer data, U.S. Geological Survey, Reston, VA, 2012.

C. Turnadge and B. D. Smerdon, “A review of methods for modelling environmental tracers in groundwater: advantages of tracer concentration simulation,” Journal of Hydrology, vol. 519, pp. 3674–3689, 2014.

B. Hagedorn, A. I. El-Kadi, A. Mair, R. B. Whittier, and K. Ha, “Estimating recharge in fractured aquifers of a temperate humid to semiarid volcanic island (Jeju, Korea) from water table fluctuations, and Cl, CFC-12 and ³H chemistry,” Journal of Hydrology, vol. 409, no. 3-4, pp. 650–662, 2011.

G. J. Houben, P. Koeniger, and J. Sültenfuß, “Freshwater lenses as archive of climate, groundwater recharge, and hydrochemical evolution: insights from depth-specific water isotope analysis and age determination on the island of Langeoog, Germany,” Water Resources Research, vol. 50, no. 10, pp. 8227–8239, 2014.

T. Röper, K. F. Kröger, H. Meyer, J. Sültenfuss, J. Greskowiak, and G. Massmann, “Groundwater ages, recharge conditions and hydrochemical evolution of a barrier island freshwater lens (Spiekeroog, Northern Germany),” Journal of Hydrology, vol. 454-455, pp. 173–186, 2012.

J. C. Vogel, “Investigation of groundwater flow with radiocarbon,” Isotopes in Hydrology, pp. 355–369, 1967.

S. Stein, A. Russak, O. Sivan et al., “Saline groundwater from coastal aquifers as a source for desalination,” Environmental Science & Technology, vol. 50, no. 4, pp. 1955–1963, 2016.

C. Andrié, P. Jean-Baptiste, C. Pierre et al., “Pore water geochemistry and mixing processes within the Tahiti barrier reef,” Geochimica et Cosmochimica Acta, vol. 62, no. 16, pp. 2809–2822, 1998.

O. Sivan, Y. Yechieli, B. Herut, and B. Lazar, “Geochemical evolution and timescale of seawater intrusion into the coastal aquifer of Israel,” Geochimica et Cosmochimica Acta, vol. 69, no. 3, pp. 579–592, 2005.

L. Bouchaou, J. L. Michelot, A. Vengosh et al., “Application of multiple isotopic and geochemical tracers for investigation of recharge, salinization, and residence time of water in the Souss–Massa aquifer, southwest of Morocco,” Journal of Hydrology, vol. 352, no. 3-4, pp. 267–287, 2008.

P. J. Stuyfzand, F. Schaars, and K. J. Van der Made, “Multitracing the origin of brackish and saline groundwaters near a dune catchment area with beach nourishment (Monster, Netherlands),” in Proceedings of the 22nd Salt Water Intrusion Meeting, pp. 93–96, Buzios, Brazil, 2012.

D. M. Han, X. F. Song, M. J. Currell, and M. Tsujimura, “Using chlorofluorocarbons (CFCs) and tritium to improve conceptual model of groundwater flow in the South Coast Aquifers of Laizhou Bay, China,” Hydrological Processes, vol. 26, no. 23, pp. 3614–3629, 2012.

E. Bryan, K. T. Meredith, A. Baker, V. E. A. Post, and M. S. Andersen, “Island groundwater resources, impacts of abstraction and a drying climate: Rottnest Island, Western Australia,” Journal of Hydrology, vol. 542, pp. 704–718, 2016.

P.-E. Oms, P. Bailly du Bois, F. Dumas et al., “Inventory and distribution of tritium in the oceans in 2016,” Science of the Total Environment, vol. 656, pp. 1289–1303, 2019.

H. Nies, I. Goroncy, J. Hermann et al., Kartierung von Tc-99, I-129 und I-127 im oberflächenwasser der Nordsee [Mapping of Tc-99, I-129 and I-127 in the surface water of the North Sea], Bundesamt für Strahlenschutz, Salzgitter, Germany, 2010.

V. E. A. Post, A. Vandenbohede, A. D. Werner, Maimun, and M. D. Teubner, “Groundwater ages in coastal aquifers,” Advances in Water Resources, vol. 57, pp. 1–11, 2013.

S. Ueda, H. Kakiuchi, K. Kondo, and J. Inaba, “Tritium concentration in fresh, brackish and sea-water samples in Rokkasho-Village, Japan, bordered by nuclear fuel cycle facilities,” Journal of Radioanalytical and Nuclear Chemistry, vol. 267, no. 1, pp. 29–33, 2005.

J. Barckhausen, “Entstehung und entwicklung der insel Langeoog [origin and development of the island of Langeoog],” Oldenburger Jahrbuch, vol. 69, pp. 239–281, 1969.

W. Voß, B. Siemon, J. Pielawa, M. Bock, and M. Ibs-von Seht, Technischer bericht hubschraubergeophysik befliegung Langeoog 2008/09 (technical report helicopter geopyhical survey Langeoog 2008/2009), Bundesanstalt für Geowissenschaften und Rohstoffe, Hannover, Germany, 2015.

S. Costabel, B. Siemon, G. Houben, and T. Günther, “Geophysical investigation of a freshwater lens on the island of Langeoog, Germany – insights from combined HEM, TEM and MRS data,” Journal of Applied Geophysics, vol. 136, pp. 231–245, 2017.

F. Bungenstock and A. Schäfer, “The Holocene relative sea-level curve for the tidal basin of the barrier island Langeoog, German Bight, Southern North Sea,” Global and Planetary Change, vol. 66, no. 1-2, pp. 34–51, 2009.

G. J. Houben, “300 years of coastal salinization research in Germany - the Homann (1718) map of the Christmas Flood of 1717,” in 25th Salt Water Intrusion Meeting (SWIM 2018), pp. 102–104, Gdańsk, Poland, 2018.

J. Sültenfuß, W. Roether, and M. Rhein, “The Bremen mass spectrometric facility for the measurement of helium isotopes, neon, and tritium in water,” Isotopes in Environmental and Health Studies, vol. 45, no. 2, pp. 83–95, 2009.

S. Lee, M. Currell, and D. I. Cendón, “Marine water from mid-Holocene sea level highstand trapped in a coastal aquifer: evidence from groundwater isotopes, and environmental significance,” Science of the Total Environment, vol. 544, pp. 995–1007, 2016.

K. R. Rushton, “Differing positions of saline interfaces in aquifers and observation boreholes,” Journal of Hydrology, vol. 48, no. 1-2, pp. 185–189, 1980.

J. H. Tellam, J. W. Lloyd, and M. Walters, “The morphology of a saline groundwater body: its investigation, description and possible explanation,” Journal of Hydrology, vol. 83, no. 1-2, pp. 1–21, 1986.

G. Marggraf, Geogene und anthropogene einflüsse auf ein isoliertes system (Düneninseln) in bezug auf die hydrochemie des grundwassers [Geogenic and anthropogenic influences on an isolated system (Dune islands) with respect to groundwater hydrochemistry], Thesis, Technische Universität Braunschweig, Braunschweig, Germany, 2005.

C. D. Langevin, D. T. Thorne Jr., A. M. Dausman, M. C. Sukop, and W. Guo, “SEAWAT version 4: a computer program for simulation of multi-species solute and heat transport,” Tech. Rep., Geological Survey (U.S.), 2008.

V. Bedekar, E. D. Morway, C. D. Langevin, and M. J. Tonkin, MT3D-USGS version 1: a U.S. Geological Survey release of MT3DMS updated with new and expanded transport capabilities for use with MODFLOW, U.S. Geological Survey, Reston, VA, 2016.

L. Stoeckl, Dokumentation der numerischen grundwasserströmungs- und transportmodellierung von Langeoog [explanatory note for the numerical groundwater flow and transportmodel of Langeoog], Federal Institute for Geosciences and Natural Resources (BGR), Unpublished report, 2018.

G. Marggraf and K. Naumann, Hydrogeologie und grundwasserdynamik (hydrogeology and groundwater dynamics). Umweltverträgliche Grundwasserbewirtschaftung in hydrogeologisch und ökologisch sensiblen Bereichen der Nordseeküste (Sustainable exploitation of groundwater in hydrogeologcally and ecologically sensitive areas along the North Sea coast), J. Petersen, R. Pott, P. Janiesch, and J. Wolff, Eds., Husum Druck- und Verlaggesellschaft, Husum, Germany, 2003.

V. E. A. Post and G. J. Houben, “Density-driven vertical transport of saltwater through the freshwater lens on the island of Baltrum (Germany) following the 1962 storm flood,” Journal of Hydrology, vol. 551, pp. 689–702, 2017.

IAEA/WMO, “Global Network of Isotopes in Precipitation. The GNIP Database,” 2016, http://www.iaea.org/water.

S. C. Doney and W. J. Jenkins, “The effect of boundary conditions on tracer estimates of thermocline ventilation rates,” Journal of Marine Research, vol. 46, no. 4, pp. 947–965, 1988.

W. Weiss and W. Roether, “The rates of tritium input to the world oceans,” Earth and Planetary Science Letters, vol. 49, no. 2, pp. 435–446, 1980.

E. H. De Jong and O. C. Swertz, Radioactieve stoffen in de zoute wateren (radioactive substances in the saline waters), Rijksinstituut voor Kust en Zee/RIKZ, Den Haag, The Netherlands, 2000.

C. Wedekind, “Tritium distribution and spreading in the North Sea and the Baltic Sea in 1980/81, as well as in the surface water of the North Atlantic in 1979,” Deutsche Hydrographische Zeitschrift, vol. 35, no. 4, pp. 177–186, 1982.

P. J. Stuyfzand, Hydrochemistry and Hydrology of the Coastal Dune Area of the Western Netherlands, KIWA N.V., Research & Consultancy Division, 1993.

A. Vandenbohede and L. Lebbe, “Effects of tides on a sloping shore: groundwater dynamics and propagation of the tidal wave,” Hydrogeology Journal, vol. 15, no. 4, pp. 645–658, 2007.

P. Pauw, P. G. B. de Louw, and G. H. P. O. Essink, “Groundwater salinisation in the Wadden Sea area of the Netherlands: quantifying the effects of climate change, sea-level rise and anthropogenic interferences,” Netherlands Journal of Geosciences - Geologie en Mijnbouw, vol. 91, no. 3, pp. 373–383, 2014.

E. M. LaBolle, G. E. Fogg, and J. B. Eweis, “Diffusive fractionation of ³H and ³He in groundwater and its impact on groundwater age estimates,” Water Resources Research, vol. 42, no. 7, article W07202, 2006.

M. R. Underwood, F. L. Peterson, and C. I. Voss, “Groundwater lens dynamics of Atoll Islands,” Water Resources Research, vol. 28, no. 11, pp. 2889–2902, 1992.

J. E. Griggs and F. L. Peterson, “Ground-water flow dynamics and development strategies at the atoll scale,” Groundwater, vol. 31, no. 2, pp. 209–220, 1993.

F. Ghassemi, J. W. Molson, A. Falkland, and K. Alam, “Three-dimensional simulation of the Home Island freshwater lens: preliminary results,” Environmental Modelling & Software, vol. 14, no. 2-3, pp. 181–190, 1998.

V. E. A. Post, A. L. Bosserelle, S. C. Galvis, P. J. Sinclair, and A. D. Werner, “On the resilience of small-island freshwater lenses: evidence of the long-term impacts of groundwater abstraction on Bonriki Island, Kiribati,” Journal of Hydrology, vol. 564, pp. 133–148, 2018.

D. K. Solomon and E. A. Sudicky, “Tritium and helium 3 isotope ratios for direct estimation of spatial variations in groundwater recharge,” Water Resources Research, vol. 27, no. 9, pp. 2309–2319, 1991.

P. M. Barlow, Ground water in freshwater-saltwater environments of the Atlantic Coast, U.S. Geological Survey, 2003.

H. H. Cooper Jr., “A hypothesis concerning the dynamic balance of fresh water and salt water in a coastal aquifer,” Journal of Geophysical Research, vol. 64, no. 4, pp. 461–467, 1959.

A. J. Smith, “Mixed convection and density-dependent seawater circulation in coastal aquifers,” Water Resources Research, vol. 40, no. 8, p. 40, 2004.

L. Stoeckl and G. Houben, “Flow dynamics and age stratification of freshwater lenses: experiments and modeling,” Journal of Hydrology, vol. 458-459, pp. 9–15, 2012.

C. I. Voss and W. W. Wood, Synthesis of geochemical, isotopic and groundwater modeling analysis to explain regional flow in a coastal aquifer of southern Oahu, Hawaii. Mathematical models and their applications to isotope studies in groundwater hydrology, IEAE-TECDOC-777, International Atomic Energy Agency (IAEA), Vienna, Austria, 1994.

Y. Yechieli, R. Yokochi, M. Zilberbrand et al., “Recent seawater intrusion into deep aquifer determined by the radioactive noble-gas isotopes ⁸¹Kr and ³⁹Ar,” Earth and Planetary Science Letters, vol. 507, pp. 21–29, 2019.

H. A. Michael, A. E. Mulligan, and C. F. Harvey, “Seasonal oscillations in water exchange between aquifers and the coastal ocean,” Nature, vol. 436, no. 7054, pp. 1145–1148, 2005.

H. Kooi and J. Groen, “Offshore continuation of coastal groundwater systems; predictions using sharp-interface approximations and variable-density flow modelling,” Journal of Hydrology, vol. 246, no. 1-4, pp. 19–35, 2001.

G. J. Houben, L. Stoeckl, K. E. Mariner, and A. S. Choudhury, “The influence of heterogeneity on coastal groundwater flow - physical and numerical modeling of fringing reefs, dykes and structured conductivity fields,” Advances in Water Resources, vol. 113, pp. 155–166, 2018.

D. Jakovovic, A. D. Werner, and C. T. Simmons, “Numerical modelling of saltwater up-coning: comparison with experimental laboratory observations,” Journal of Hydrology, vol. 402, no. 3-4, pp. 261–273, 2011.

C. Lu, P. K. Kitanidis, and J. Luo, “Effects of kinetic mass transfer and transient flow conditions on widening mixing zones in coastal aquifers,” Water Resources Research, vol. 45, no. 12, 2009.

L. M. Bexfield, B. C. Jurgens, D. M. Crilley, and S. C. Christenson, Hydrogeology, water chemistry, and transport processes in the zone of contribution of a public-supply well in Albuquerque, New Mexico, 2007-9, U.S. Geological Survey, Reston, VA, 2012.