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

Authors: Vincent E. A. Post; Georg J. Houben; Leonard Stoeckl; Jürgen Sültenfuß
Journal: Geofluids
 (2019)

---

## Abstract

Tritium (3H) and its daughter product 3He 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 3H and 3He 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 3H concentration of the North Sea. It is found that most bomb-related tritiogenic 3He still resides in the freshwater lens, making it a useful tracer for young (<60 years) groundwater. Differences in dispersive transport between 3H and 3He can cause an apparent age bias on the order of 10 years. Under favourable conditions, 3H 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 3H in the transition zone does not originate from seawater but from freshwater affected by the bomb peak. Yet in models with a low (αL=0.5m) dispersivity, no 3H was sequestered into the transition zone and the transition zone width was underestimated. Better results were obtained with αL=5m, 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.


---

## Body

## 1. Introduction

Tritium (3H) 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 (3He) is simultaneously analysed and after the so-called tritiogenic contribution (3Hetrit) is separated from other 3He 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 3Hetrit over the 3H concentration and the half-life constant:
(1)τ=1λln1+H3etritH3,in which λ=0.05626/yr is the radioactive decay rate of 3H. 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:
(2)R=nBτlnBB−d,where n is the porosity, B (L) is the lens thickness, and d (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 3H 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 3H in brackish groundwater to the mixing of young groundwater with old, 3H-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 3H levels in the North Sea [18]. As a result, Stuyfzand et al. [14] were able to attribute 3H values of up to ~30 TU (Tritium Units) in saline groundwater to the intrusion of 3H-rich seawater. The main motivation for the present work was the discovery of 3H 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 3H-enriched seawater, similar to the findings by Stuyfzand et al. [14].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.
(a)
(b)
(c)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 3H within the freshwater but also on the fate of 3H 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 3H in the transition zone (i.e., whether it originates from the freshwater or the seawater). In what follows, the notation 3He is used to indicate the tritiogenic 3He 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 (Figure1). 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 m3/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 m3/yr (Figure 2), while by 2011 the pumped volume was down to 333,000 m3/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 (Q) from the freshwater lens (blue bars) as well as seasonally distributed abstraction per unit of surface area (Q′) versus time (red line); (b) calculated groundwater recharge (R, bars) and 3H concentration of groundwater recharge (3Hrech, solid red line) and seawater (3Hsea, 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 y axis in (b) indicates the mean recharge rate (0.87mm/d=319mm/yr) 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; 1TU=3H/1H ratio of 10−18=0.118TU) was achieved by using the 3He accumulation method. For all sampling sites, duplicate samples were taken and analysed.Transition zone water samples were obtained from observation well 17 (Figure1), which has a long screen (from 33.6 to 57.6 m below sea level) reaching into the brackish water. Profiles of specific conductance (σ25, 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)3H concentration, and (d) 3He 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 σ25 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 αL=5m. 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 (Figure4). SEAWAT [31] was used to solve the following mass balance equations for fluid and solute, respectively,
(3)∂nρ∂t+∇⋅ρq→=ρ¯qs,(4)∂nC∂t=∇⋅nD→⋅∇C−∇⋅q→C−qsCs+Γ,where ρ is the fluid density (M L-3), q→ is the Darcy flux (L T-1), ρ¯ is the density of water entering (or leaving) through a source (or sink) [M L-3], qs is the source’s volumetric flow rate per unit volume of aquifer (T-1), C is the concentration (M L-3), D→ is the hydrodynamic dispersion tensor (L2 T-1), and Cs is the solute concentration associated with water sources and sinks (M L-3). The hydrodynamic dispersion tensor D→ is a function of the longitudinal (αL) and transverse (αT) dispersivity (L), as well as the effective molecular diffusion coefficient (De) (L2 T-1).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 (Q) and recharge (R) time series used are shown in Figures 2(a) and 2(b), respectively. The input time series for 3H in recharge (3Hrech) and seawater (3Hsea) 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-3 T-1) in equation (4) accounts for the effects of radioactive decay. Three dissolved species, Cl, 3H, and 3He, were simulated. For Cl, Γ=0, while for 3H, the decay rate is Γ=−λCH, where the subscript H is used to indicate that C represents the 3H concentration. For 3He, Γ=λCH, i.e., the growth rate of tritiogenic 3He is equal in magnitude but opposite in sign to the decay rate of 3H. 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 (Table1). 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 (A31)
A3
A4
A5
A6
A7
A8
A9
Transient recharge
X
X
X
X
X
X
X
X
AnisotropicK
X
X
X
X
X
X
X
LayeredK
X
X
X
X
X
X
Clay layer
X
X
Discontinuous clay layer
X
X
Abstraction
X
X
XThe effect of the vertical anisotropy on the3H 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 2<Kx/Kz<20. This number was used as a subscript to indicate the simulation, e.g., A35 is for simulation A3 with Kx/Kz=5. Since simulation A2 differs only from A3 in that it has Kx/Kz=1, it has been equivalently labelled as A31. Some simulations with a constant 3H 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 rainfall3H 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 3H 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.The3H 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 3H 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 the3H 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 3H=0.2TU 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 of1.2×106m2. 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 αL=0.3m, 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 αL=0.5m (with αL/αT=10). This value is intermediate between the values used by Stuyfzand [41] and Vandenbohede and Lebbe [42], who adopted values of αL=0.1m and αL=0.2m, respectively, for the freshwater lenses in the Dutch and Belgian coastal dunes, and the value used by Pauw et al. [43] who used αL=2m 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 5×10−10, 6.7×10−10, and 2×10−9m2/s were used for Cl, 3H, and 3He, respectively [44].


## 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; 1TU=3H/1H ratio of 10−18=0.118TU) was achieved by using the 3He accumulation method. For all sampling sites, duplicate samples were taken and analysed.Transition zone water samples were obtained from observation well 17 (Figure1), which has a long screen (from 33.6 to 57.6 m below sea level) reaching into the brackish water. Profiles of specific conductance (σ25, 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)3H concentration, and (d) 3He 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 σ25 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 αL=5m. 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 (Figure4). SEAWAT [31] was used to solve the following mass balance equations for fluid and solute, respectively,
(3)∂nρ∂t+∇⋅ρq→=ρ¯qs,(4)∂nC∂t=∇⋅nD→⋅∇C−∇⋅q→C−qsCs+Γ,where ρ is the fluid density (M L-3), q→ is the Darcy flux (L T-1), ρ¯ is the density of water entering (or leaving) through a source (or sink) [M L-3], qs is the source’s volumetric flow rate per unit volume of aquifer (T-1), C is the concentration (M L-3), D→ is the hydrodynamic dispersion tensor (L2 T-1), and Cs is the solute concentration associated with water sources and sinks (M L-3). The hydrodynamic dispersion tensor D→ is a function of the longitudinal (αL) and transverse (αT) dispersivity (L), as well as the effective molecular diffusion coefficient (De) (L2 T-1).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 (Q) and recharge (R) time series used are shown in Figures 2(a) and 2(b), respectively. The input time series for 3H in recharge (3Hrech) and seawater (3Hsea) 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-3 T-1) in equation (4) accounts for the effects of radioactive decay. Three dissolved species, Cl, 3H, and 3He, were simulated. For Cl, Γ=0, while for 3H, the decay rate is Γ=−λCH, where the subscript H is used to indicate that C represents the 3H concentration. For 3He, Γ=λCH, i.e., the growth rate of tritiogenic 3He is equal in magnitude but opposite in sign to the decay rate of 3H. 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 (Table1). 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 (A31)
A3
A4
A5
A6
A7
A8
A9
Transient recharge
X
X
X
X
X
X
X
X
AnisotropicK
X
X
X
X
X
X
X
LayeredK
X
X
X
X
X
X
Clay layer
X
X
Discontinuous clay layer
X
X
Abstraction
X
X
XThe effect of the vertical anisotropy on the3H 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 2<Kx/Kz<20. This number was used as a subscript to indicate the simulation, e.g., A35 is for simulation A3 with Kx/Kz=5. Since simulation A2 differs only from A3 in that it has Kx/Kz=1, it has been equivalently labelled as A31. Some simulations with a constant 3H 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 rainfall3H 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 3H 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.The3H 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 3H 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 the3H 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 3H=0.2TU 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 of1.2×106m2. 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 αL=0.3m, 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 αL=0.5m (with αL/αT=10). This value is intermediate between the values used by Stuyfzand [41] and Vandenbohede and Lebbe [42], who adopted values of αL=0.1m and αL=0.2m, respectively, for the freshwater lenses in the Dutch and Belgian coastal dunes, and the value used by Pauw et al. [43] who used αL=2m 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 5×10−10, 6.7×10−10, and 2×10−9m2/s were used for Cl, 3H, and 3He, respectively [44].


## 4. Results and Discussion

### 4.1. Transient Effects

The results of simulations A1 and A2 are shown in Figure4. 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 the3H bomb peak can clearly be recognized as a quasihorizontal, elongated zone of high 3H 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 3H concentrations were at their highest, groundwater recharge rates were relatively low (Figure 2(b)). Therefore, the total mass of 3H that entered the freshwater lens was considerably lower and resulted in a lower total 3H mass value for simulation A2 than for simulation A1 (Figure 6). For example, at the peak in September 1966, total 3H contained within the freshwater lens amounted to 1.4×108Bq for simulation A2, compared to 1.8×108Bq for simulation A1.Figure 5
Cross sections showing contour plots of the simulated3H 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 (=A31), and (c) A320. The dashed white line in (a) and (b) indicates the position of the vertical profiles shown in Figure 7.
(a)
(b)
(c)Figure 6
Temporal evolution of the total3H in Bq (solid lines) and 3He (dashed lines) in the freshwater lens (groundwater with less than 2.5% seawater) for simulations A1 (red) and A2 (blue).The graph in Figure6 further shows the total 3He produced from tritium in the freshwater lens. The maximum was reached during the mid-1980s but, as 3He production in the freshwater lens fell because of the declining 3H concentrations, the 3He lost by discharging groundwater started to exceed production and total 3He has been decreasing since then. Nevertheless, most (about 80%) of this 3He was still present within the freshwater lens by the end of the simulation. Groundwater dating back to the 1960s thus remains identifiable by high 3He concentrations it acquired due to the 3H peak in recharge at that time. This provides opportunities to use 3He in lieu of 3H 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 3H in rainfall being variable in time on the apparent age τ. Tritium enters the saturated zone with the recharge, but the decay product 3He forms in situ in the aquifer. When the input signal of 3H is steady, steady-state 3He and 3H 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 3H is variable in time, such as in the case of the 3H bomb peak, 3He 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 3He concentrations. As the concentration gradients of both 3H and 3He 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 3H 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 3H 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 3He concentration.Figure 7
(a) Apparent ageτ (solid lines) and chloride concentration (dashed lines) for simulations A1 and A2 and (b) the 3H (solid lines) and 3He (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 3H input of 9 TU (as indicated by the subscript c) and for the 3H 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 3He/3H ages even when the 3H input signal is constant in time. In the case of simulation A2c, the maximum difference of the 3He/3H 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 3H input.The apparent age shift due to the nonconstant3H input is also apparent when the sum of 3H and 3He is plotted as a function of τ (equation (1)), which is often done to identify mixing effects in field data (e.g., [8]). If 3H and its daughter 3He were to remain within the same parcel of groundwater, their sum would reflect the 3H concentration at the time of recharge. In that case, a groundwater sample would have a matching data point somewhere on the atmospheric 3H 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 3H and 3He and its effect on τ, the peak of the summed modelled 3H and 3He concentrations of the groundwater does not coincide with the peak of the atmospheric 3H 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 the3H concentration in precipitation with time (same data as shown in Figure 2(b)) to the sum of the simulated and measured 3H and 3He concentrations in groundwater versus τ for simulations A1 (red dots), A2 (blue dots), and A7 with αL=5m (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 elevated3H 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 3H-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 3H 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 3H (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 Kx/Kz be 1, 2, 5, 10, and 20 by varying Kz. 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 Kx and the recharge R, 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 ~600<x<~800m.In contrast to the relative constancy of the freshwater lens geometry, Figure5(c) clearly shows that the decrease of Kz significantly reduces the maximum penetration depth of the funnel-shaped body of high 3H saltwater (measured by zmin, the elevation of the deepest point where 3H>0.02 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 Kx/Kz changed from 1 to 20. The effect on the inflow of seawater across the seafloor (Qs) is shown in Figure 9. The shape of the relationship between Qs and Kz is the same as the relationship between Qs/Qf and V∗ defined by Smith [52] (his Figure 2), where V∗=Kz/qf×ρs−ρf/ρf. In the setup used by Smith [52], Qf and qf 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 Kz changed between the simulations and all other parameters that determine V∗ remained constant.Figure 9
Dependency of (a) seawater inflowQs and (b) zmin (the elevation of the deepest point where 3H> 0.02 TU) on vertical hydraulic conductivity Kz for simulation A3 with different anisotropy ratios. The corresponding anisotropy ratios are displayed along the top horizontal axis of the graph. Qs 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 Qs and zmin, 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 use3H 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, 3H measured in water samples from a multilevel observation well at the shoreline might be used to infer the magnitude of Qs. A relatively high 3H 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 (Figure3), 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 3H 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.Figure 10
Contour plots of the simulated3H 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)
(b)
(c)
(d)
(e)
(f)A complex distribution of chloride and3H 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 3H 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 and3H 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 Figures10(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 3H 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, 3H activities were generally lower in the simulations with abstraction. This is attributed partially to the removal of water containing 3H 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 the3H and 3He 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 3H 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 3H 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 3He 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 (Figure4) 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 (Figure10(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 αL=5m (i.e., ten times the original value, and αT 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 3H and 3He 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 did3H originating from intruded seawater reach the location of the observation well, which most likely rules it out as a source of the observed 3H 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 αL=5m did result in elevated 3H concentration in the transition zone, suggesting that the mechanism put forward by Sivan et al. [12], i.e., freshwater containing bomb-derived 3H mixes with 3H-free seawater, can be responsible for 3H in the transition zone in coastal aquifers.But even with a longitudinal dispersivityαL=5m, the Langeoog models underestimated the measured 3H 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 αL=5m were not considered but may be required to explain the width of the transition zone.


## 4.1. Transient Effects

The results of simulations A1 and A2 are shown in Figure4. 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 the3H bomb peak can clearly be recognized as a quasihorizontal, elongated zone of high 3H 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 3H concentrations were at their highest, groundwater recharge rates were relatively low (Figure 2(b)). Therefore, the total mass of 3H that entered the freshwater lens was considerably lower and resulted in a lower total 3H mass value for simulation A2 than for simulation A1 (Figure 6). For example, at the peak in September 1966, total 3H contained within the freshwater lens amounted to 1.4×108Bq for simulation A2, compared to 1.8×108Bq for simulation A1.Figure 5
Cross sections showing contour plots of the simulated3H 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 (=A31), and (c) A320. The dashed white line in (a) and (b) indicates the position of the vertical profiles shown in Figure 7.
(a)
(b)
(c)Figure 6
Temporal evolution of the total3H in Bq (solid lines) and 3He (dashed lines) in the freshwater lens (groundwater with less than 2.5% seawater) for simulations A1 (red) and A2 (blue).The graph in Figure6 further shows the total 3He produced from tritium in the freshwater lens. The maximum was reached during the mid-1980s but, as 3He production in the freshwater lens fell because of the declining 3H concentrations, the 3He lost by discharging groundwater started to exceed production and total 3He has been decreasing since then. Nevertheless, most (about 80%) of this 3He was still present within the freshwater lens by the end of the simulation. Groundwater dating back to the 1960s thus remains identifiable by high 3He concentrations it acquired due to the 3H peak in recharge at that time. This provides opportunities to use 3He in lieu of 3H 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 3H in rainfall being variable in time on the apparent age τ. Tritium enters the saturated zone with the recharge, but the decay product 3He forms in situ in the aquifer. When the input signal of 3H is steady, steady-state 3He and 3H 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 3H is variable in time, such as in the case of the 3H bomb peak, 3He 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 3He concentrations. As the concentration gradients of both 3H and 3He 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 3H 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 3H 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 3He concentration.Figure 7
(a) Apparent ageτ (solid lines) and chloride concentration (dashed lines) for simulations A1 and A2 and (b) the 3H (solid lines) and 3He (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 3H input of 9 TU (as indicated by the subscript c) and for the 3H 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 3He/3H ages even when the 3H input signal is constant in time. In the case of simulation A2c, the maximum difference of the 3He/3H 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 3H input.The apparent age shift due to the nonconstant3H input is also apparent when the sum of 3H and 3He is plotted as a function of τ (equation (1)), which is often done to identify mixing effects in field data (e.g., [8]). If 3H and its daughter 3He were to remain within the same parcel of groundwater, their sum would reflect the 3H concentration at the time of recharge. In that case, a groundwater sample would have a matching data point somewhere on the atmospheric 3H 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 3H and 3He and its effect on τ, the peak of the summed modelled 3H and 3He concentrations of the groundwater does not coincide with the peak of the atmospheric 3H 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 the3H concentration in precipitation with time (same data as shown in Figure 2(b)) to the sum of the simulated and measured 3H and 3He concentrations in groundwater versus τ for simulations A1 (red dots), A2 (blue dots), and A7 with αL=5m (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 elevated3H 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 3H-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 3H 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 3H (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 Kx/Kz be 1, 2, 5, 10, and 20 by varying Kz. 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 Kx and the recharge R, 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 ~600<x<~800m.In contrast to the relative constancy of the freshwater lens geometry, Figure5(c) clearly shows that the decrease of Kz significantly reduces the maximum penetration depth of the funnel-shaped body of high 3H saltwater (measured by zmin, the elevation of the deepest point where 3H>0.02 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 Kx/Kz changed from 1 to 20. The effect on the inflow of seawater across the seafloor (Qs) is shown in Figure 9. The shape of the relationship between Qs and Kz is the same as the relationship between Qs/Qf and V∗ defined by Smith [52] (his Figure 2), where V∗=Kz/qf×ρs−ρf/ρf. In the setup used by Smith [52], Qf and qf 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 Kz changed between the simulations and all other parameters that determine V∗ remained constant.Figure 9
Dependency of (a) seawater inflowQs and (b) zmin (the elevation of the deepest point where 3H> 0.02 TU) on vertical hydraulic conductivity Kz for simulation A3 with different anisotropy ratios. The corresponding anisotropy ratios are displayed along the top horizontal axis of the graph. Qs 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 Qs and zmin, 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 use3H 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, 3H measured in water samples from a multilevel observation well at the shoreline might be used to infer the magnitude of Qs. A relatively high 3H 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 (Figure3), 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 3H 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.Figure 10
Contour plots of the simulated3H 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)
(b)
(c)
(d)
(e)
(f)A complex distribution of chloride and3H 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 3H 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 and3H 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 Figures10(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 3H 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, 3H activities were generally lower in the simulations with abstraction. This is attributed partially to the removal of water containing 3H 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 the3H and 3He 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 3H 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 3H 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 3He 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 (Figure4) 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 (Figure10(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 αL=5m (i.e., ten times the original value, and αT 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 3H and 3He 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 did3H originating from intruded seawater reach the location of the observation well, which most likely rules it out as a source of the observed 3H 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 αL=5m did result in elevated 3H concentration in the transition zone, suggesting that the mechanism put forward by Sivan et al. [12], i.e., freshwater containing bomb-derived 3H mixes with 3H-free seawater, can be responsible for 3H in the transition zone in coastal aquifers.But even with a longitudinal dispersivityαL=5m, the Langeoog models underestimated the measured 3H 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 αL=5m 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 3H 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 3H and its daughter 3He in the aquifer.The model results can be synthesised in terms of groundwater age classes as shown in Table2, 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 3He was found to still be present within the freshwater lens, it may be possible to use 3He as a tracer to identify freshwater recharged around the time of the bomb peak (a role traditionally fulfilled by 3H). Because 3He 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 groundwater3H and 3He concentration (all numbers in TU) for simulation A7 with αL=5m.
3H<0.2
0.2<3H<123He<10
0.2<3H<123He≥10
3H>12
Freshwater
≥60 yr (3He ≈ 5)
<10 yr
≥20 yr
<10 yr (3He<5)
Saltwater
≥40 yr (3He ≈ 0.2)
Mixture of old seawater andseawater<20yr
Mixture of seawater and ≥20 yr freshwater (brackish)
<10 yr (3He<12)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 simulated3H 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 αL=5m 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 the3H 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 3H 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 elevated3H 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 3H can be found in the subsea portion of the aquifer, as well as on submarine groundwater discharge patterns. Therefore, groundwater 3H 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 3H or 3He variations can be resolved with sufficient analytical precision.



---

*Source: 1494326-2019-10-07.xml*