Introduction

Barrier islands are highly dynamic and unique environments, which provide many ecosystems services and form attractive tourist destinations. Fresh groundwater often is the only water resource. It is stored as a thin veneer in freshwater lenses surrounded by more saline water, making groundwater abstraction a delicate trade-off between meeting water demand and preventing salinization. The freshwater lens condition and the amount of water that can be safely withdrawn from it depend strongly on groundwater recharge. Knowing the recharge rate is thus of utmost importance to ensure sustainable use.

The exposure of barrier islands to sea currents and wave actions makes their coastlines susceptible to erosion. This, and inundations by storm floods, puts the fresh groundwater at direct risk (Post and Houben 2017). The erosion threat often prompts local authorities to take coastal defence measures. Typical measures include the heightening of existing dunes or even creating artificial coastal dunes by spraying and bulldozing dredged sand in the form of elongated bodies along the foredunes (Caljé et al. 2018; Thorenz 2008).

While the creation of artificial dunes and islands is becoming more widespread (Huizer et al. 2016; van Ginkel 2015), little is known about groundwater recharge processes in these systems. The initial scarcity of vegetation limits transpiration and allows a greater fraction of precipitation to reach the water table and recharge the freshwater lens. Over the course of several years, the artificial dunes undergo a change from an almost unvegetated “white dune” type to a “brown dune” type, vegetated by mosses, grass, and low shrubs. Studies in the coastal dunes in the Netherlands have shown that under bare or poorly vegetated dunes, recharge occurs during the whole year, with the fraction of rainfall becoming recharge being lowest in May and June and highest in November and December (Stuyfzand 1993). Under densely vegetated dunes, recharge only occurs from August/September to April, with a soil-moisture deficit developing during the remaining months.

The objective of the present work is to determine groundwater recharge rates and temporal dynamics for two contrasting sites on Langeoog Island in northern Germany: an old, vegetated dune and a recently created, sparsely vegetated dune. To this aim, the study involved soil-water stable isotope (2H, 18O) measurements that were simulated using a new version of the HYDRUS-1D code recently extended by Zhou et al. (2021) to include isotope fractionation. Soil-moisture isotope measurements have been used successfully to constrain recharge rates in many studies (Koeniger et al. 2016; Tewolde et al. 2019; Mattei et al. 2020), and models of varying complexity have been applied for this purpose (Barbecot et al. 2018; Boumaiza et al. 2021). Recharge estimates can also be based on chloride concentration measurements, but given that chloride inputs are highly variable in space and time close to the coast, e.g., due to sea spray, isotopes form a valuable alternative, albeit with added complexity due to measurement procedures and fractionation effects. The results from this work form a demonstration example of the application of the new HYDRUS-1D code for recharge estimation more generally. They also provide insight into the effects of artificial dune creation on groundwater recharge to freshwater lenses in coastal areas.

Study area

The island of Langeoog is an elongated sandy barrier island located at the northern rim of the shallow intertidal Wadden Sea, off the coast of northern Germany (Fig. 1). Undifferentiated Holocene and Pleistocene (glacio-)fluvial sediments form the base of the island, overlain by Holocene tidal flat deposits of the Wadden Sea (at depths of 10–20 m below sea level) and, finally, Holocene beach and dune sands (Bungenstock and Schäfer 2009). The island formed around 2800–2200 years ago, with the first dunes forming at around AD 100. Due to constant aeolian sand movement, vegetation was very scarce, and settlement only began in the thirteenth century, after some soil had formed (Barckhausen 1969). The dunes grew considerably between the thirteenth and seventeenth centuries. Catastrophic storm floods, especially the Christmas flood of 1717, separated the dune belt into three individual parts (Houben 2018). Since the beginning of the twentieth century, the dunes were repeatedly widened and heightened with sand dredged from the North Sea and stabilized after deposition by planting beach grass (marram).

Fig. 1
figure 1

a–b Maps showing the location of Langeoog Island and the soil sampling sites and weather stations; c–d soil surface at profile sites LB01 (older dune) and LB02 (younger dune), respectively. Photograph in d courtesy of M. Howahr, OOWV. Base map and data from OpenStreetMap and OpenStreetMap Foundation (OpenStreetMap 2021)

Full size image

Due to generally favourable conditions in the dune area (sandy soils, low transpiration and interception losses by the predominantly low vegetation), three freshwater lenses formed below the island. Groundwater extraction in the western lens began in 1909 with three wells located in the western Kaap dunes, close to the village of Langeoog (Fig. 1). They were abandoned in 1989. The main well field was started in 1938 in the Heerenhus dunes, more to the east of the village but became fully operational in 1957. While the local population is small (2,000 people), intensive tourism is the main driver of water demand, using up to three-quarters of the total extraction of around 350,000 m3 year−1 (Houben et al. 2014). Demand is highly seasonal, with consumption during the main tourist season up to eight times higher than in the off-season.

The German barrier islands have been intensively studied because of the importance of their freshwater lenses for water supply. Detailed descriptions of their hydrogeology, including Langeoog, can be found in Costabel et al. 2017; Holt et al. 2021; Houben et al. 2014; Post et al. 2019; Post and Houben 2017; Röper et al. 2012. Several previous studies have attempted to quantify groundwater recharge on Langeoog Island. Seliger (2002) reported groundwater recharge estimates based on three different calculation procedures. Two of these were based on a calculation of actual evapotranspiration, either on a daily or annual basis. The third method calculated the percolation rate using a soil-water balance model, based on inputs of monthly rainfall and potential evapotranspiration. Depending on the method, the 13-year average value calculated by Seliger (2002) ranged from 300 to 563 mm year−1. For individual years, the differences in calculated recharge rates between the methods amounted up to a factor of 3.

Houben et al. (2014) used groundwater ages based on tritium-helium data to calculate recharge rates and found that rates varied across the island depending on land use and vegetation. The highest recharge rates of 400 mm year−1 or more were found below the dune valleys (Houben et al. 2014). Age-depth profiles from wells below dune crests yielded recharge rates lower by a factor of two. In the built-up village and the forested part of Langeoog, groundwater recharge rates were found to be as low as 100 mm year−1, the former caused by soil sealing, the latter by leaf interception and transpiration.

Liebscher (1970) published the results of lysimeter studies conducted on the nearby island of Norderney, which is located some 20 km west to Langeoog, between 1956 and 1967. Percolation rates for an unvegetated lysimeter in the dune area covered with grass ranged between 346 mm (1959, total rainfall of 511 mm) and 939 mm (1966, total rainfall of 1,070 mm). Overall, the reported annual percolation rate was 86% of the annual rainfall sum (Liebscher 1970).

Data and methods

Soil-water samples

Soil samples were collected from two sites adjacent to the Pirolatal on the northern edge of the dune belt (Fig. 1) on 25 June 2019. Profile LB01 (old dune) is located on the flank of a natural dune vegetated with grasses and herbaceous vegetation (Fig. 1d) at an elevation of 8.46 m above mean sea level (amsl) as determined by differential global positioning system (GPS) measurements. Profile LB02 (young dune) is located on top of an artificial white dune (elevation 9.25 m amsl) with scant vegetation (Fig. 1e). The sands there were deposited for the first time in 2008 and again in 2015.

The soil samples were obtained by hand augering using a Dutch auger system, and collected in laminated coffee bags (WEBER Packaging, Güglingen, Germany) that were sealed using a Ziplock system as well as by heat-sealing. They were transported the next day to the laboratory, where they were stored refrigerated until extraction. Chloride concentrations was analyzed by ion chromatography (Dionex ICS-3000) after aqueous extraction from oven-dried (105 °C for 24 h) soil samples according to the procedures in the EN 15934:2012–11 standard. For stable isotope analysis the soil water was extracted using a modified cryogenic vacuum extraction method (Koeniger et al. 2011). For each soil sample, three subsamples were taken for chemical and isotope analysis.

Stable isotope analyses were conducted at the stable isotope laboratory of BGR in Hannover (Germany) using a Picarro L2120-i water vapour analyser connected to a vaporizer and autosampler system. Isotope values are given in the standard delta notation per mil (‰) versus VSMOW. The data sets were corrected for machine drift during the run and normalized to the VSMOW/SLAP scale by assigning a value of 0‰ and –55.5‰ (δ18O) and 0‰ and –428‰ (δ2H) to VSMOW and SLAP, respectively. External analytical reproducibility for the laser spectrometer, defined as the standard deviation of a control standard during all runs, was better than 0.5 and 1.5‰ after soil-water extraction (2-sigma standard deviation) for δ 18O and δ2H, respectively (Koeniger et al. 2011).

Soil-moisture contents were determined gravimetrically by oven-drying at 105 °C as well as by weighing the original and cryogenically dried sample weights. The gravimetric soil-water content was converted to a volumetric soil-water content using a dry bulk density of 1.6 g cm−3 for quartz sand with a porosity of 0.4. Grain size analyses were conducted using a Camsizer laser particle sizer (Retsch, Germany) after sieving the samples to remove the fraction <63 μm (constituting only 0.2% by weight on average for all samples) and >2 mm (absent in all samples).

Soil hydraulic properties and wettability

Due to the dry conditions, undisturbed soil core samples could not be taken at the sites of the soil-water profiles because the sand was too loose. Instead, soil hydraulic properties were obtained from soil cores collected at two sites at the nearby waterworks (Fig. 1), sampled and analysed during an earlier project in 2014. Water characteristic curve and unsaturated hydraulic conductivity were measured in the BGR soil laboratory in Hannover with the evaporation method (Wind 1969) using a Hyprop device (Meter Group, Munich, Germany). Eight 250-cm3 samples, collected at 6–11 cm depth, were taken to the laboratory and saturated with tap water before the evaporation measurement (630 and 82 data pairs, respectively). The saturated hydraulic conductivity was determined using a laboratory permeameter (12 replications; Eijkelkamp, Giesbeek, The Netherlands) on 250-cm3 soil cores sampled at 50–55-cm depth.

The wettability of the dune sand was measured with the water-drop-infiltration-time test (WDPT). This test can be performed in situ or on disturbed samples in the laboratory. Dekker and Ritsema (1994) introduced the concept of actual and potential water repellency. Measurements on moist/wet soil material indicate actual water repellency (WDPTact), i.e., water repellency of the soil material at a certain water content (typically determined in situ). As the wettability of soil material is dependent on its water content, the WDPT was also determined on oven-dried (40 °C) samples to measure the potential water repellency (WDPTpot). Potential water repellency is assumed to be the maximum water repellency of a sample (even though a higher drying temperature might create higher WDPTpot values in some samples).

In situ measurements were carried out along a transect of about 100 m, crossing a dune valley and the corresponding dune crest (Fig. 1a). Three drops of water (approx. 50 μL) were placed on a smoothed, bare (unvegetated) spot of dune sand, and the time it took to infiltrate was recorded (WDPTact). At each location, approximately 100 g of dune sand was sampled, which was dried at 40 °C for 7 days in the laboratory. After equilibration to room temperature, three drops of water were placed on the smoothed sample surface, and the time was recorded for the drop to infiltrate into the sand (WDPTpot). Measured WDPT-values were classified according to Bisdom et al. (1993).

Numerical modelling

The processes in the unsaturated zone at sites LB01 and LB02 were simulated using a modified version of the HYDRUS-1D software package that is able to simulate isotope fractionation by evaporation (Zhou et al. 2021). The model considers water flow with snow hydrology, evaporation, root water uptake, and heat and solute transport. Variably saturated water flow was described using the Richards equation:

$$ \frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left[K(h)\frac{\partial h}{\partial z}+K(h)\right]-S $$
(1)

where θ is the volumetric water content, z is the depth below the soil surface (L), t is time, K(h) is the hydraulic conductivity (L T−1), which is a function of the pressure head (h) (L), and S (T−1) is a sink term representing root water uptake.

The water retention characteristic curve θ(h) and the hydraulic conductivity function K(h) are described as suggested by van Genuchten (1980) and Mualem (1976):

$$ \theta (h)=\left\{\begin{array}{c}{\theta}_{\mathrm{r}}+\frac{\theta_s-{\theta}_{\mathrm{r}}}{{\left(1+{\left|\alpha h\right|}^n\right)}^m}\\ {}{\theta}_{\mathrm{s}}\end{array}\right.\kern1.75em {\displaystyle \begin{array}{c}h<0\\ {}h\ge 0\end{array}} $$
(2)
$$ K(h)={K}_{\mathrm{s}}\frac{{\left[1-{\left(\alpha h\right)}^{n-1}{\left[1+{\left(\alpha h\right)}^n\right]}^{-m}\right]}^2}{{\left[1+{\left(\alpha h\right)}^n\right]}^{\frac{m}{2}}} $$
(3)

where θs and θr are the saturated and residual water contents, respectively, α (L−1), n, and m (= 1 – 1/n) are empirical parameters, and Ks is the saturated hydraulic conductivity (LT−1). These parameters were determined by fitting Eqs. (2) and (3) to the measured soil hydraulic property data.

The site 1 data were therefore used to parameterise the van Genuchten-Mualem model in HYDRUS-1D. An uncertainty analysis was conducted to explore the sensitivity of the yearly recharge rates to the parameter uncertainty. This involved running 100 realisations of each model (i.e., LB01 or LB02) in which θs, θr, α and Ks were randomly sampled from a normal distribution (see Table S1 in the electronic supplementary material (ESM1).

Because winter precipitation is sometimes in the form of snow, snow cover was modelled using the simplified approach implemented in HYDRUS-1D (Šimůnek et al. 2008). This approach assumes that precipitation is in the form of snow when air temperatures are below zero °C and, consequently, does not infiltrate immediately into the soil during freezing periods. Snowmelt is simulated to be proportional to air temperature above freezing, using a snow melting constant of 0.43 cm day−1 K−1 (i.e. the HYDRUS-1D default).

Root water uptake was calculated using the Feddes et al. (1978) water-uptake-reduction model using the parameters provided by Voortman (2018) for a similar dune vegetation site in the Netherlands, which has comparable climatic and soil conditions as the present study site. For profile LB01 (natural dune, vegetated), the leaf area index (LAI) was set to 1, which is at the lower range of values for stands of dune grass and herbaceous vegetation reported of a similar dune area in Belgium (Willaert et al. 2005). The LAI for profile LB02 (artificial dune, scant vegetation) was treated in a special manner discussed later in the article. Based on observations made during augering, roots were assumed to be present in the upper 40 cm of the soil profiles. The root distribution function in HYDRUS-1D [b′(x)] was assumed to decrease exponentially with depth to reflect the fact that root density is highest just below the land surface.

For water flow simulations, an “atmospheric boundary condition with surface runoff” and a “free drainage” boundary condition (BC) were used for the upper and lower boundaries, respectively. The upper BC involves alternatively Neumann and Dirichlet boundary conditions, depending on conditions in the soil. Within a prescribed pressure head interval, the Neumann BC is used. Once the limit of this interval is reached, the Dirichlet BC is used. The “free drainage” BC is a zero-pressure-head-gradient-boundary condition. No immobile water regions were considered in the simulations.

HYDRUS-1D solves the advection-dispersion equation for solute transport. The dispersivity was optimised to fit the measured isotope values. Cauchy and Neumann BCs were used for solute transport at the top and bottom boundaries of the soil profile, respectively. Different to the standard HYDRUS-1D code, the Cauchy BC at the top was modified to consider fractionation by evaporation, as described by Zhou et al. (2021). The isotopic composition (δ2H, δ18O) of the evaporation flux is calculated according to the Craig-Gordon model (δE, CG) as

$$ {\delta}_{\mathrm{E},\mathrm{CG}}=\left[{\mathrm{RH}}_{\mathrm{s}}\bullet {\alpha}_{\mathrm{i}}^{\ast}\bullet {\delta}_{\mathrm{i}}^{\mathrm{ls}}-{\mathrm{RH}}_{\mathrm{a}}^{\prime}\bullet {\delta}_{\mathrm{a}}-{\mathrm{RH}}_{\mathrm{s}}\bullet {\varepsilon}^{\ast }-{\varepsilon}_k\right]/\left({\mathrm{RH}}_{\mathrm{s}}-{\mathrm{RH}}_{\mathrm{a}}^{\prime }+{\varepsilon}_{\mathrm{k}}/{10}^3\right) $$
(4)

where \( {\delta}_{\mathrm{i}}^{\mathrm{ls}} \) and δa (‰) are the isotopic composition of the soil water (ls = liquid/soil), and the atmospheric vapour, respectively, RHs is the relative humidity at the soil surface [−], \( {\mathrm{RH}}_{\mathrm{a}}^{\prime } \) is the relative humidity of the atmosphere [−] normzalized to the soil surface temperature, \( {\alpha}_{\mathrm{i}}^{\ast } \) and \( {\alpha}_{\mathrm{i}}^{\mathrm{total}} \) are the equilibrium fractionation factor [−] and the total fractionation factor [−], respectively, and ε∗ and εk are the equilibrium and kinetic fractionation enrichment (‰), respectively. Details about the calculation procedure and model parameters can be found in Zhou et al. (2021). Isotope fractionation is implemented for the first (upper) model cell. Root water uptake does not cause fractionation; hence, any water that is lost by plant transpiration does not lead to isotopic enrichment of the soil moisture.

Weather data

Daily precipitation data measured by the Deutscher Wetterdienst (DWD) from Langeoog (station number 02861) for the period between 1941 and 1994 were downloaded from the climate data centre (DWD 2021). Daily precipitation, wind velocity, air humidity, mean air temperature, minimum air temperature, maximum air temperature, and sunshine duration are available from 1947 onward for the barrier island Norderney (station number 3631). The DWD data also include the potential evapotranspiration (ETo) according to Penman-Monteith from 1991 onward, which were also used in the numerical model. The island’s water supply company (Oldenburgisch-Ostfriesischer Wasserverband, OOWV) operates a weather station in the dune area on Langeoog (see Fig. 1 for location). The measured daily precipitation and temperature were used as input in the numerical models, with missing data substituted by the DWD data from Norderney.

Precipitation samples for isotope analysis were collected at the OOWV weather station (Fig. 1) as biweekly to monthly cumulated samples using a precipitation isotope sampler (PALMEX RS1) (Gröning et al. 2012) from July 2012 and June 2019. This period was chosen as the HYDRUS-1D model simulation period. The delta values of the atmospheric water vapour (δ2Ha, δ18Oa) were determined by assuming they are in equilibrium with the prevailing rainfall values.

Lysimeter data

The former municipal water supply company of Langeoog operated a lysimeter on the island of Langeoog from July 1963 to October 1970 and from November 1974 to April 1977. In total, data for 118 months are available, including the summer of 1976, which was particularly hot and dry. The archived data for the lysimeter were supplied by the current water supplier (OOWV) in paper form. The unvegetated lysimeter had a surface area of 1 m2, and the inferred depth was 1.45 m. It was located in the dunes near the old well field, situated west of the village (Fig. 1). The soil had a porosity ranging between 0.39 and 0.43, a coefficient of uniformity (d60/d10) in a narrow range between 1.1 and 1.3, and a hydraulic conductivity range of 2.9 ∙ 10−5–4.0 ∙ 110−5 m s−1. These parameters indicate a well-sorted sediment, typical for dune sands. The method of determining the hydraulic conductivity was not documented.

Monthly sums of precipitation data were also reported as part of the lysimeter dataset; however, these showed significant differences compared with the monthly sums of the DWD data. There is a good correlation between the monthly sums for the first series of lysimeter data (July 1963 to October 1970), with the lysimeter values generally being higher by 15%. This may reflect spatial variability of the rainfall or a difference in instruments used. For the second series of lysimeter data (November 1974 to April 1977), there is a very poor correlation, suggesting a problem with the reported monthly sums for the lysimeter. This part of the data series was therefore not considered. ETo data (required for modelling purposes) prior to 1991 were not available, nor were wind speed data required to calculate ETo available for the period July 1963 to October 1970. Instead, the vapour pressure deficit (VPD) was determined from temperature and relative humidity data, which were available for the entire measurement period. Linear regression was then made between ETo and VPD for the 1991 and beyond data, and this regression relationship was used to calculate ETo based on VPD for July 1963 to October 1970.

Results

Figures 2 and 3 show the measurements for both soil profiles from June 2019. These will be discussed in detail in the following subsections.

Fig. 2
figure 2

Depth distribution for profile LB01 of a the volumetric water content (θ), b the stable isotope delta values of the soil moisture, c the chloride concentration of the soil moisture, and d the median grain size

Full size image
Fig. 3
figure 3

Depth distribution for profile LB02 of a the volumetric water content (θ), b the stable isotope delta values of the soil moisture, c the chloride concentration of the soil moisture, and (d the median grain size. The sample between 30 and 40 cm depth has a chloride concentration of 1,110 mg L−1 and is not visible

Full size image

Soil-moisture content, isotope values, and chloride concentrations

Soil-moisture contents show only minor variability and are close to the residual water content θr (discussed in more detail further on). The lowest values were measured in the upper 40 cm of profile LB02. There is a good correspondence between the oven-dried and cryogenically extracted values, providing confidence that the cryogenic extraction method provided representative water samples for the isotope analysis. The large error bars in Fig. 2a are due to the measurement of triplicate samples, and the high variability indicates that the isotope measurements of these samples potentially have to be ignored because of incomplete extraction.

The isotope delta values of profile LB01 show only little variability below a depth of 1 m. Above this level, the samples are more enriched in the heavier isotopes. The variations are somewhat larger in the profile of profile LB02 (especially for δ2H), and strong enrichment of the topmost sample is evident from the high delta values. These samples also have very high chloride values. The sample between the 30 and 40 cm depth of profile LB02 had a chloride concentration of 1,110 mg L−1 and is not shown in the graph. A clear chloride concentration gradient can be seen from the soil surface down to a depth of 80 cm in profile LB01.

The scatter plot in Fig. 4a reinforces the correlation between δ18O and chloride that is visible in the soil profiles, especially for profile LB01. In the scatter plot of δ2H versus δ18O (Fig. 4b), most of the data points straddle along the local meteoric water line (LWML, δ2H = 7.23 × δ18O + 4.56, R2 = 0.97), except for the most enriched samples, which are located towards the right of the LMWL, which is evidence for evaporation.

Fig. 4
figure 4

a δ18O versus the chloride concentration and b δ2H versus δ18O of the two soil-moisture profiles. The rainfall data points included in b were analysed in samples collected at the OOWV weather station from July 2012 and June 2019

Full size image

Soil physical parameters

The median grain size of both profiles varies within a narrow range, which is generally between 170 and 190 μm for LB01 and between 190 and 210 μm for LB02 (Figs. 2 and 3). For all samples, the grain size fractions 112–200 μm and 200–355 μm combined make up 90% by weight of the total, indicating a uniform sediment. The particle size distribution curves for all soil samples are shown in Fig. S1 in ESM1.

The van Genuchten-Mualem model parameters (Eqs. 2 and 3) were determined by an unconstrained inverse parameter estimation to the 630 data points from the Hyprop measurements; however, the resulting value of n was higher than 5, and this caused frequent crashes of the HYDRUS-1D model. Trial runs demonstrated that values of n = 4.8 resulted in stable model runs, so parameter estimation was subsequently conducted using this value. The resulting parameters are shown in Fig. 5. The two sets of data that were fitted were the samples at site 1 and the surface samples at site 2, located near the Langeoog waterworks. All samples come from the Heerenhus dunes and are thus of similar origin and lithology. The site 2 samples at the 50–55 cm depth interval are very similar to the surface samples and were not separately fitted.

Fig. 5
figure 5

a Soil-water retention curves and b relationship between K(h) and the log soil matric potential pF. The blue curve was fitted to the data points of site 1, the orange curve to the surface sample at site 2. The corresponding van Genuchten-Mualem parameter values are provided next to the curves. The green-shaded areas indicate the range of retention and K(h)/pF relationships considered in the uncertainty analysis

Full size image

The measured saturated hydraulic conductivity varied by more than one order of magnitude (between 1.8 ∙ 10−5 and 2.2 ∙ 10−4 m/s, or 156 and 1,963 cm/day), with the lower conductivity values occurring in the uppermost part of the unsaturated zone. Only the four samples that were analysed using a tension-infiltrometer yielded values that were in a pF range that allowed constraining the K(h) relation. The resulting value of l = 4,030 is extremely high compared to the recommended value of 0.5 (Mualem (1976)). Although this has been reported for other datasets as well (Yates et al. 1992), the value was not used in the subsequent HYDRUS-1D simulations because it resulted in extremely long runtimes. Moreover, trial simulations showed that the parameter set based on the site 2 data resulted in much poorer model fits to the isotope data than the site 1 data.

The dune sand on Langeoog Island shows severe to extreme water repellency when dry, as indicated by high WDPTpot values (Fig. 6). The in-situ measurement showed lower WDPTact times compared to WDPTpot times. Water repellency was not at its maximum at the time of measurement; however, actual water repellency was observed at the dune crest. This is not only due to potentially lower water contents of the sand at these locations, as WDPTpot times also increased at these locations. A slight trend can be seen that the dune crest is more water repellent compared to the dune valley.

Fig. 6
figure 6

Potential and actual water repellency as determined with the WDPT test on samples from a transect across a sand dune (see the location in Fig. 1). Repellency classes according to Bisdom et al. (1993)

Full size image

Model results

The simulated isotope profiles are compared to the measurements in Fig. 7b–e. The fit between the model that considers fractionation and the data is especially good in the deeper part of the soil profiles for both sites LB01 and LB02. Because the van Genuchten-Mualem parameters were determined independently from the soil hydraulic property measurements (Fig. 5), the only parameter used to obtain the fit in Fig. 7 was the dispersivity. The adjusted value was 4 cm for LB01 and 2 cm for LB02, which is within the reported range of other studies (2–20 cm) with comparable transport distances (80–200 cm; Stumpp et al. 2012; Vanderborght and Vereecken 2007).

Fig. 7
figure 7

a Rainfall (P), potential evapotranspiration (ET), and rainfall δ2H values versus time as input to the model (δ18O values not shown). Measured and simulated δ2H values for b profile LB01 and c profile LB02. Measured and simulated δ18O values for d profile LB01 and e profile LB02

Full size image

The animated version of Fig. 7 (see ESM2) shows the seasonal dynamics of the isotope profiles. The animation starts 250 days before the soil samples were taken at a moment when the delta values are highest in the top part of the profile due to evaporation during the summer of 2018. The water with these high delta values gets pushed down as rainfall increases during the winter. The animation clearly shows how the downward migration accelerates during periods of high recharge. While the bulk of the enriched soil moisture exits through the bottom model boundary during the winter months of 2019, slightly elevated delta values remain visible in the model as well as the data at a depth just below 4 m. The peak that is the result of precipitation in January 2019 with high delta values remains preserved, albeit attenuated, at a depth of 3 m. As the top part of the profile starts drying out from May 2019 onward, the delta values there start to increase by evaporative enrichment, eventually leading to the profile as observed during the fieldwork. The high values persist near the land surface as long as there is no net downward flux to flush out the enriched soil moisture.

For both models, the upper parts of the measured isotope data could not be fitted without considering fractionation in both models (Fig. 7). This was expected because the samples were taken in summer, and the preceding months had been unusually dry and warm. Interestingly, the model without fractionation did not fit the data points in the deeper part of profile LB01, but it did so for profile LB02. No matter what parameter combinations were attempted, models with LAI = 1 systematically overpredicted the measured delta values in the deeper part of the profile at LB02 if fractionation was considered. The animated version of Fig. 7 makes it clear that the bulges and troughs in the measured delta values between 2 and 4.5 m depths formed during October 2018 and March 2019, i.e., during the winter months. From this, it was inferred that direct evaporation during this period is apparently negligible at LB02. In order to mimic this behaviour, the model was run with LAI = 10−6 between April and September and LAI = 100 from October to March. These parameter values have no physical meaning, but were merely chosen to partition potential evapotranspiration into transpiration (so that fractionation is effectively switched off) in autumn and winter and into evaporation in spring and summer. The result is shown as the solid line in Fig. 7b. The interpretation in terms of physical processes will be discussed further on in the article.

The model’s ability to match the observed isotope patterns provides confidence that the simulated water fluxes are accurate. The temporal variability of the calculated recharge rates is shown in Fig. 8. The summer (May–October) and winter (November–April) sums are listed in Table 1. Recharge rates tend to be highest in winter, but recharge occurs throughout the year, even during the driest summer months. The differences between sites LB01 and LB02 are small in winter, but the rates are higher for the less vegetated LB02 profile than for LB01 in summer. The yearly sums are included in Table 2 for comparison with other published studies.

Fig. 8
figure 8

a Monthly rainfall and potential evapotranspiration (ETo) versus time and b side by side comparison of monthly recharge rates (R) for sites LB01 and LB02

Full size image
Table 1 Six-month sums for precipitation (P), potential evapotranspiration (ETo), recharge (R), and the ratio of recharge and precipitation for sites LB01 and LB02. The values are grouped by winter and summer values
Full size table
Table 2 Yearly rainfall (P), potential evapotranspiration (ETo), recharge (R), and the ratio of recharge and precipitation for sites LB01 and LB02
Full size table

The results of the uncertainty analysis are summarised in Figs. S2 and S3 and Table S2 in ESM1. It can be seen that the standard deviations for the yearly recharge rates are small, ranging from 9 to 17 mm year−1 for LB01 to 5 to 13 mm year−1 for LB02. This demonstrates that the calculated annual recharge rates are not very sensitive to the soil physical parameters. On the other hand, the shape of the isotope profiles is sensitive to the soil physical parameters, as shown by Fig. S3 in ESM1.

Lysimeter

Just like the modelling results, the 1963–1970 lysimeter data were split into 6-month winter and summer periods. Table 3 shows that percolation amounted to between 29 and 53% of precipitation during summer, while in winter, this ranges between 85 and 93%. Similar to the HYDRUS-1D modelling results for the soil profiles, the data show that recharge is not restricted to the cooler winter months.

Table 3 Six-month sums for precipitation (P), potential evapotranspiration (ETo), recharge (R), and the ratio of recharge and precipitation for the lysimeter. The values are grouped by winter and summer values
Full size table

As a further validation of the HYDRUS-1D model, the same soil hydraulic parameters were used to calculate the reported percolation rates from the lysimeter. As the lysimeter was unvegetated, root water uptake was not considered. The results in Fig. 9 show that there is a very good agreement up until the summer of 1968. After this point in time, the modelled values become systematically higher than the reported values. Not enough information about the experimental setup is available to identify the cause for this deviation. Nevertheless, the good fit during the first 5 years provides additional confirmation that the rates calculated by the model are representative of the dune sand on the island.

Fig. 9
figure 9

a Monthly rainfall and potential evapotranspiration (ETo) versus time and b observed and modelled cumulative percolation rates for the lysimeter. The cumulative values are shown to emphasise the systematic deviation after the summer of 1968 (b)

Full size image

Figure 10 compares the summer and winter sums for the 1963–1970 lysimeter data and the 2012–2019 model results. It is clear that the data points follow similar trends. It is important to bear in mind that the lysimeter was only 1.45 m deep, while the modelled soil profiles had a length of 4.5 m. A comparison on a monthly basis, therefore, looks much less favourable because infiltrated water needs much longer to travel across the 4.5-m than across the 1.45-m soil column. The similar trends of the 6-month sums provide confidence about the accuracy of these recharge rates.

Fig. 10
figure 10

Winter (Nov–Apr) and summer (May–Oct) sums of recharge versus precipitation for the 1963–1970 lysimeter data as well as the model-based values for LB01 and LB02 for the period 2012–2019. An ‘exclamation mark’ marks the data points representing the summer of 1968 and later, i.e., the period during which the lysimeter data started to diverge from the model simulation

Full size image

Discussion

The preceding results show that the soil-moisture isotope profiles can be modelled accurately by the HYDRUS-1D code if fractionation is considered. Over the entire depth of the profiles, the delta values in profile LB01 are somewhat higher than at LB02. This is somewhat unexpected because the vegetation is more abundant at LB01 than at LB02, and hence transpiration, which does not cause fractionation, might be expected to be more important. One reason could be that the downward displacement at LB01 is slower than at LB02 due to root water uptake, especially during summer months, so that there is more opportunity for direct evaporation, and hence fractionation can impart a measurable enrichment of the soil moisture.

Interestingly, the data of profile LB02 could only be simulated with the model by assuming that the soil moisture did not experience fractionation during the winter season, but was significantly affected by direct evaporation during the summer. The explanation for this is thought to lie in the fact that during the unusually dry and warm summer of 2019, a significant part of the vegetation died off (see also Fig. 1e), thus increasing the importance of direct evaporation over transpiration. By setting the LAI to a very low value in summer, this behaviour could be replicated by the model. A comparison of aerial images from July 2015 (with 407 mm of summer rainfall, Table 1) and August 2019 clearly shows a scarcer vegetation cover at LB02 (see ESM1). Voortman et al. (2015) contended that vegetation dieback is an important process that must be considered in recharge estimates of coastal dune areas. It is not clear if this effect is equally important every year. Only a single snapshot in time of the soil-water isotope profile is currently available for this study area, which is representative for a period of about 9 months prior to sampling. Repeated measurements alongside modelling could help in testing this hypothesis and provide more accurate recharge estimates.

Vegetation dieback is likely to enhance soil-water repellency by exposing the sand surface, which causes the sand to become extremely dry. Measurements showed significant water repellency of the dune sands, especially on the dune crest (Fig. 6). The difficulty of water to infiltrate would expose it to direct evaporation, which would impart enrichment of the heavier isotopes. Whether or not this contributes to the enrichment of the soil moisture depends on whether the water on the surface evaporates completely or only partially, and if it infiltrates with the next rain event that wets up the soil sufficiently so that infiltration can occur again.

To what extent observed soil-water repellency affects the recharge processes could not be ascertained due to a lack of data. In water repellent soils, infiltration can only occur when the ponding height (and thus the water pressure) exceeds the resistance to flow due to water repellency (Bauters et al. 2000). The microtopography of the surface therefore plays an important role in this. It can lead to the formation of preferential flow with an uneven distribution of soil-moisture contents and downward flow rates. The process is favoured under initially dry conditions, and when it occurs, water can move vertically downward along the preferential flow paths at greater rates than under wet conditions. Any rainfall that is unable to infiltrate into the repellent soil may be evaporated (reducing recharge) or flow down the dune slope (redistributing recharge). Houben et al. (2014) suggested that the inferred difference in recharge rates between dune crests and valleys was due to overland flow and accumulation and subsequent infiltration of water near the dune foot.

Scanlon et al. (2002) advocated using multiple techniques to enhance the confidence of groundwater recharge estimates. While the lysimeter data and modelled soil-water profile rates showed a good resemblance in this study, the calculated recharge rates are higher than values reported in previous studies on Langeoog Island. For example, the 13-year average that was reported by Seliger (2002) was only 224–420 mm year−1, depending on the calculation method used, which corresponds to 31–58% of the mean annual precipitation. These estimates, however, were based on theoretical calculations that involved numerous assumptions and uncertain parameter values and were not constrained by any soil or groundwater data. Houben et al. (2014) used groundwater ages (based on 3H/3He measurements of groundwater samples) in combination with the Vogel (1967) model to determine the spatial variation of mean annual recharge rates. They found values between 230 and 400 mm year−1 in the dune area, with the highest values reported for the dune valleys. These values are a lot lower than the annual mean value for profile LB01 (527 mm year−1, Table 2), which is the more representative of the two sites studied here for the dune area as a whole.

The explanation for this difference might lie in the fact that, despite the presence of vegetation, the plant water use is lower at LB01 than in the bulk of the dune area due to differences in vegetation type. While shrubs are absent at profile LB01, they are abundant in other parts of the dune area. Studies in the Netherlands and Belgium have demonstrated that their presence enhances evapotranspiration (Samson et al. 2005; Stuyfzand 1993). Moreover, some plants, like marrum, are known to grow their roots all the way down to the water table and can take water directly from the saturated zone. This additional loss of water is not considered in the models of this study, which are restricted to the unsaturated zone. Alternatively, the Houben et al. (2014) rates could be an underestimate, possibly because the flow conditions in the aquifer are more complex than assumed in their application of the one-dimensional Vogel (1967) model to calculate recharge based on groundwater age. It should also be recalled that recharge rates based on lysimeter measurements and soil-water profiles represent point estimates, whereas rates based on groundwater age data are affected by processes integrated over a certain area and time span.

The lysimeter data reported by Liebscher (1970) for a dune on the island of Norderney would be expected to be in a similar range as this study’s estimates. The reported percolation rates are 75% of the annual rainfall (period 1956–1967). The same percentage was reported by Stuyfzand (1993) based on a giant (625 m2) lysimeter in the coastal dunes of the Netherlands that was filled with soil that was almost bare, except for a scant moss cover. Therefore, the percentages found, based on the HYDRUS-1D model for sites LB01 and LB02 (Table 2), could well be realistic, and the differences with the studies by Houben et al. (2014) are due to different methodologies or uncaptured processes like direct groundwater use by vegetation or overland flow. Figure 11 summarizes the results from both lysimeters, the models for both Langeoog sites and the age data.

Fig. 11
figure 11

Comparison of recharge rates obtained from the unvegetated lysimeters in Norderney (N) and Langeoog (L), the HYDRUS-1D models for sites LB01/LB02 and the age dating by Houben et al. (2014). Symbols represent individual years. Legend for bars, see LB02 model

Full size image

The model’s ability to simulate the observed soil-moisture isotope values demonstrates that this method holds a lot of promise to improve recharge estimates of coastal dune areas. While in the current study, the isotope contents were determined by destructive sampling, new technologies make a nondestructive experimental setup possible as well. Using portable laser-based measurement instruments that measure the soil-water vapour isotope content in-line, a high-resolution time series of soil-water isotope contents can be derived. High-frequency sampling of the rainfall isotope content and actual measurements of atmospheric water vapour isotopic composition could lead to improvements of the model and better insight into the temporal dynamics of the soil-water and recharge processes.

Knowledge of the dynamic nature of groundwater recharge is especially important given the unknown impact of climate change. More intense rainfall events may lead to more infiltration, but could also enhance overland flow. It may be anticipated that periods of prolonged drought and warmer summers will intensify soil-water repellency and vegetation dieback. Studies in the Netherlands (Witte et al. 2012) have shown that the proportion of moss and bare soil increases during summer droughts, increasing recharge in a warmer climate due to decreased transpiration. Intensified monitoring in combination with modelling is thus needed to establish how recharge processes will respond to climate change in the future.

Conclusions

This study investigated the recharge rates and processes below two contrasting dune types on the Langeoog barrier island in northern Germany. The recharge is crucial for sustaining the freshwater lens that is exploited for water supply and determines the hydrology, and thereby the ecology, of the dune area. The two investigated sites differ in their genesis and vegetation type. The dune of profile LB01 is an older dune that is covered by grass and herbaceous vegetation, whereas the dune of profile LB02 is a recent dune that was created to reinforce the island’s coastal protection and is scarcely vegetated.

Soil samples were collected from the unsaturated zone by hand auguring down to the water table. The grain size and soil-moisture content of the sand were analysed. Soil-water samples for stable water isotope analyses were obtained by cryogenic extraction, and chloride concentrations of the soil water were also measured. The results show that evaporation induces fractionation of the isotopes in the upper part of the profile, which was expected, as the samples were taken during the warm and dry summer of 2019.

Recharge rates were calculated using a new version of HYDRUS-1D, modified to simulate isotope fractionation during evaporation, and fitting the simulated isotope-depth profiles to the measurements. The model outcomes confirm the importance of fractionation, and it was found that fractionation effects are more pronounced at the profile with the most vegetation (LB01). This is attributed to slower downward percolation of the soil moisture as water is lost to root water uptake, providing more opportunity for direct evaporation in the upper part of the soil profile. For profile LB02, the model could be matched to the data by assuming that transpiration dominates in winter versus direct evaporation in summer. The inferred absence of transpiration in summer is tentatively attributed to vegetation dieback, but a dataset spanning multiple years is needed to test this assertion.

The calculated recharge rates are consistent with lysimeter measurements that were made on Langeoog Island in the past. Both the model results and the lysimeter data of this study show that recharge occurs year-round and that, as expected, the winter months contribute most of the recharge. The recharge rates are significantly higher than the rates reported by Houben et al. (2014) that were based on groundwater age data. The discrepancy is attributed to potential methodological issues and, predominantly, the role of vegetation. The presence of dune shrub is likely to result in lower recharge rates than those found at sites LB01 and LB02. Further research is needed to establish the significance of other processes, like soil-water repellency and overland flow.

It is expected that the methodology applied in this study can yield important insights into the dynamical nature of the recharge processes when measurements are taken repeatedly and frequently. Such studies are not only needed at the present study site, but across a range of geographical locations as recharge processes will respond to climate change in presently unknown ways.