The W-E trending barrier island of Langeoog is located at the German North Sea coast, along the northern rim of the intertidal Wadden Sea (Fig. 1a). The water supply on the island relies exclusively on local groundwater resources, present as three freshwater lenses located underneath sand dune bodies. The western lens is the only one currently in use for water supply. A detailed description of the hydrogeology and water resources management of Langeoog is given by Houben et al. (2014). Freshwater lenses are fragile bodies and prone to over-exploitation, storm floods and sea-level rise. Numerical groundwater models are useful tools to predict and counteract the effects of such events. The quality of its results, however, depends strongly on the accuracy of the input data, e.g. information on the geological structures and hydraulic parameters in the subsurface. A limited number of boreholes provides only sparse input data considering the complexity of the aquifer systems. The aim of this study is to show how airborne and ground-based geophysical methods can supplement geological information.
In substantial parts of Northern Germany the Late Elsterian Lauenburg clay acts as aquiclude. On Langeoog, however, it was affected by later glacial push processes (Streif, 1990, Streif, 2004), which commonly render it hydraulically ineffective. Therefore, the permeable Pleistocene and Holocene layers generally act as one aquifer. On a local scale, however, especially in the vicinity of the production wells, the presence of local clay layers is expected to have a significant impact on the movement of groundwater. As a consequence, knowledge about the local spatial distribution of aquitards is required for a sustainable groundwater management.
Geophysical methods have the potential to identify local aquicludes and aquitards. On Langeoog, several ground-based geophysical surveys using the vertical resistivity sounding (VES) method were conducted in the years 1970, 1976, 1983, and 1995, comprehensively summarized by Tronicke (1997). An interesting finding from these studies is an increase of the freshwater lens thickness northeast of the managed area. Tronicke (1997) assumed that this is the result of the precautions implemented in the 1980s, e.g. the artificial raising of dunes to protect the lens from overwash. Another possible reason is a decrease of groundwater extraction after 1980 (Houben et al., 2014). However, the VES results did not provide indications of clay layers in the area of investigation. In 2008 and 2009, helicopter-borne electromagnetic (HEM) surveys were performed at the North Sea coast of Germany including the island of Langeoog (Siemon et al., 2015). This dataset is used in this study as the basis to delineate the current dimensions of the freshwater resources on the island. However, the results of resistivity methods (HEM, TEM, VES) are, by their nature, ambiguous regarding lithology and groundwater salinity (e.g. Vouillamoz et al., 2007, Klimke et al., 2013, Siemon et al., 2015).
In addition to HEM, a ground-based geophysical campaign in 2012 and 2013 was conducted in cooperation with the local water supplier Oldenburgisch-Ostfriesischer Wasserverband (OOWV). The main objective of these investigations was to supplement the lithological information available from boreholes in order to improve the hydrogeological basis for the numerical modelling of the freshwater lens in the western part of the island. The survey consisted mainly of transient electromagnetic (TEM) and magnetic resonance sounding (MRS) measurements. This combination of methods has been proven to be successful in similar studies: Vouillamoz et al. (2012) used combined interpretation of TEM and MRS data for retrieving hydraulic properties of a thin freshwater lens on a barrier island. Behroozmand et al. (2012) jointly inverted these data using common layer models. Similarly, Günther and Müller-Petke (2012) used VES instead of TEM on the neighbouring island of Borkum, which exhibits a very similar geological setting.
Down to the depths considered in this study (< 60 m), the island of Langeoog comprises three principal geological units (Barckhausen, 1969): the base is formed by undifferentiated glaziofluvial sediments, mostly of Pleistocene age. They contain the Lauenburg clay, formed at the end of the Elsterian as a meltwater basin sediment. Its depth (15 to 35 m below sea level) and thickness (typically a few meters) can vary considerably. It is absent in parts of the island and this indicates a strong influence by later glacial push processes, probably during the Drenthe stage of the following Saalian glaciation. The Pleistocene is overlain by Holocene Wadden Sea deposits which can be found down to depths of 10 to 20 m below sea level. The youngest deposits are dune and beach sands found at and above the current sea level (Bungenstock and Schäfer, 2009). The island itself is geologically very young. About 2800 to 2200 years ago it formed as a sand bar, onto which later dune sands accumulated (Barckhausen, 1969). The island remained largely free of vegetation for several centuries due to constant Aeolian sand movement (Barckhausen, 1969). It was not until the 13th century before enough soil had formed to allow human settlement. Until the 17th century the dunes considerably grew by Aeolian sand accumulation. Several catastrophic storm surges, especially around the 1720s, caused large breaches of the dune belt and flooding of the fresh groundwater bodies. In the following centuries, vegetation of the dunes, man-made earthworks and coastal protection measures established the three dune complexes, which today dominate the morphology of the island: west (settled area), middle (Melkhörn) and east.
Fig. 1b shows the area of investigation and Fig. 2 shows two geological cross sections based on borehole information (A-A′ and B-B′) that demonstrate the general hydrogeological setting of Langeoog. The coarse sand at the bottom of the profiles corresponds to the Pleistocene, overlain by a clayey to silty aquitard layer (Lauenburg clay) with intercalated sand layers. The aquifer in use comprises the fine to medium grained Holocene sands and the fine dune sands. Fig. 2b shows that the aquitard can be absent locally.
The positions of the HEM flight lines and the TEM/MRS measurements are depicted in Fig. 1b together with the locations of observation boreholes and production wells. Detailed information on boreholes and geophysical measurement are given in Table 1. While the HEM survey covers the entire island, the ground-based geophysical methods focused on the eastern part of the managed lens. This was because the borehole information in this area is very sparse, whereas the borehole density in the middle and in the western part is sufficient regarding a comprehensive geological interpretation.
Helicopter-borne electromagnetics (HEM)
HEM is an airborne electromagnetic method which uses systems that are towed below a helicopter. The frequency-domain HEM system (Resolve, Fugro Airborne Surveys) is part of the airborne geophysical system of the Federal Institute for Geosciences and Natural Resources (Germany) consisting of electromagnetics, magnetics and radiometrics, which was applied successfully in many groundwater exploration surveys (e.g. Siemon et al., 2007, Siemon et al., 2009). The HEM device, with both horizontal coplanar (HCP) and vertical coaxial (VCX) transmitter-receiver coil configurations, is operated at six frequencies ranging from 0.4 to 133 kHz. The signals measured are the relative secondary magnetic fields of eddy currents induced in a conductive ground by the transmitter coils. At each frequency, they are individually picked up by the corresponding receiver coils, sampled every 0.1 s and split into their in-phase (I) and out-of-phase (Q) components, relative to the transmitter signal. Standard HEM data processing is applied comprising calibration, drift correction, levelling and correction for artificial electromagnetic (EM) disturbances (Siemon et al., 2009). The latter is particularly necessary for regions where man-made effects, caused by electric installations or urban areas, distorted HEM measurements which then have to be erased. The resulting gaps can be closed by calculating synthetic data values from interpolated half-space parameter grids (Siemon et al., 2011).
The airborne survey over the island of Langeoog was part of a larger survey in Eastern Friesland conducted in 2008 and 2009 (Siemon et al., 2015). The nominal separation of the N-S lines was 250 m and the W-E tie lines were 2000 m apart. The distance between consecutive values was in the order of 3.5 m assuming an average flight speed of 125 km/h during the survey. The height of the sensor system was kept 30 to 40 m above ground surface wherever possible.
The HEM data (only HCP data which are suitable for investigating layered structures) were inverted to layered-earth models using an iterative Marquardt-Levenberg inversion process (Sengpiel and Siemon, 2000, Siemon et al., 2009). The starting models necessary for the inversion were derived from apparent resistivity vs. centroid depth sounding curves (Siemon, 2001, Siemon, 2012) for both block and smooth inversion with 6 and 19 layers, respectively. While the block inversion calculated all model parameters (here: 6 layer resistivities and 5 layer thicknesses), the 18 layer thicknesses below the cover layer were fixed for the smooth inversion. In order to avoid oscillations of the 19 layer resistivities, the factor controlling the damping of the model update for the smooth inversion was chosen about three times stronger than for the block inversion. The layer thicknesses increasing with depth were individually derived from the sounding curves at each position enabling variable model depths with respect to the penetration of the HEM fields. The maximum depth of investigation, which depends on the local conductivity of the subsurface, was about 60 m below sea level on the island of Langeoog due to the existence of saltwater below the freshwater lenses limiting the further penetration of the HEM fields.
The HEM inversion results were displayed as resistivity maps at specified depths and vertical resistivity sections showing the 1D resistivity models along a survey profile including the topographic relief in m above mean sea level (asl).
Transient electromagnetics (TEM)
Time-domain or transient electromagnetics (TEM) is a technique commonly used in mining and groundwater exploration (Nabighian and Macnae, 1991). Usually, transmitter loops of squared shape with a side length of 25 to 150 m are laid on the ground (e.g. Fitterman and Stewart, 1986, Auken et al., 2003). A direct current is applied and cut-off abruptly in the transmitter circuit. The cut-off process initiates a change in the magnetic field of the transmitter coil, which causes the propagation of eddy currents into the subsurface. The current flow in the subsurface produces a secondary magnetic field that is measured using a receiver coil at the surface. The shape of the transient signal is a measure of the subsurface resistivity distribution. This is reconstructed by an appropriate inversion technique. If shallow structures up to a depth of a few decameters are to be resolved, the induced voltage in the receiver loop must be measured as early as possible after shutting down the transmitter, i.e. a measurement system with a short cut-off time of less than 10 μs for the transmitting current is required (Barsukov et al., 2007).
For the TEM measurements we used a rectangular single loop configuration and loop sizes of 25 and 50 m side length. Because the target depth was less than 60 m, we used a TEM device (TEM-Fast) optimized for shallow investigations with a cut-off time of 3 μs (Barsukov et al., 2007). A 50-Hz notch filter was applied to suppress anthropogenic EM noise. Inconsistent data points, i.e. outliers, were erased before inversion. As for HEM, smooth and block inversion schemes were applied. The TEM smooth inversion was conducted using the same 19 fixed depth layers as for the HEM data for a comparison of both methods. The starting model was a homogeneous halfspace with a resistivity of 100 Ωm. The smoothness constraints for the TEM inversion were realized by including a derivation matrix (de Groot-Hedlin and Constable, 1990) into the least squares (LS) scheme so that the roughness of the resistivity distribution is minimized along with the data misfit. Different levels of smoothness are realized by a factor weighting the two terms in the inversion process. The regularization factor was manually determined such that the final inversion result represents the smoothest model that is able to fit the measurement data without bias. For the block inversion, no smoothness constraints were defined but a Marquard-Levenberg damping as used in the HEM data inversion. The parameters of the starting model for the block inversion, i.e. the number of layers (between three and five) as well as resistivities and thicknesses of the layers, were determined after analysing the general trend of the subsurface resistivity distribution given by the smooth inversion result. Resistivity and thickness of each layer are independent parameters. After the data misfit has reached a LS minimum, a bootstrapping algorithm was applied to estimate the uncertainty of each parameter, i.e. each parameter was subsequently varied in both directions until the misfit between the original and varied model response exceeded the data misfit (Günther and Müller-Petke, 2012).
Magnetic resonance sounding (MRS)
Magnetic resonance sounding (MRS) is a ground-based hydrogeophysical method based on the principle of nuclear magnetic resonance (NMR). A transmitter loop, similar in shape and geometry to the loop used for TEM, is placed on the surface and is energized by an alternating current matching the Larmor frequency of the spins of the protons in the water molecules (e.g. Legchenko and Valla, 2002, Yaramanci and Müller-Petke, 2009, Behroozmand et al., 2014). The corresponding electromagnetic field propagates into the ground and excites the proton spins, which are forced to flip away from their equilibrium positions in the Earth’s magnetic field. After shutting down the transmitting current, the spins relax back, which results in a measureable signal observed by the receiver loop as a function of time. The MRS signal is an exponentially decaying (relaxing) signal that provides information on the volumetric water content (θMRS) from the initial amplitude (linear proportionality) and the pore geometry from the relaxation behaviour (e.g. Kenyon, 1997, Mohnke and Yaramanci, 2008, Grunewald and Knight, 2011, Dlugosch et al., 2013). An increasing pore size yields an increasing relaxation time (T2⁎). Thus, estimates of hydraulic conductivity can be obtained from MRS, if adequate calibration data is available, i.e. T2⁎- θMRS -permeability relationship from pumping tests or laboratory investigations (Legchenko et al., 2004, Mohnke and Yaramanci, 2008). The depth range of an MRS measurement is increased by increasing the pulse moment, which is the product of applied current and duration of excitation (= pulse length). Normally, an MRS data set consists of an ensemble of single measurements with varying pulse moments covering a certain depth range. The sensitivity function of this ensemble is calculated considering the loop layout and the subsurface resistivity distribution (Legchenko and Valla, 2002, Behroozmand et al., 2014), which is normally provided by resistivity methods such as VES (Günther and Müller-Petke, 2012) or TEM (Behroozmand et al., 2012). The vertical distributions of θMRS and T2⁎ are reconstructed by inversion (e.g. Mohnke and Yaramanci, 2008, Mueller-Petke and Yaramanci, 2010).
Two limitations need to be noted: First, water in small pores, e.g. in clayey and silty sediments, might already be relaxed after the dead time of the MRS measurement and does not participate in the recorded signal in this case. Therefore, aquitards are often indicated by reduced or vanishing θMRS rather than small T2⁎. The effective dead time considers the instrumental dead time, additional time delays given by filtering procedures during post-processing (e.g. bandpass-filter), and also relaxation effects during the excitation (Walbrecker et al., 2009, Dlugosch et al., 2011). Second, the site-specific calibration necessary for a quantitative estimation of hydraulic conductivity demands hydraulic data from pumping tests or laboratory investigations on aquifer material or, better still, drill cores. If these are not available, calibration must be based on literature values, which can cause unknown and unidentifiable systematic errors. On the other hand, θMRS can be estimated without any calibration (e.g. Legchenko and Valla, 2002, Legchenko et al., 2004, Yaramanci and Müller-Petke, 2009).
The hydraulic conductivity K is usually computed from MRS by an empirical equation:
where a, b, c are empirical parameters to be calibrated. However, most commonly b = 2, while for a either 1 or 4 are suggested (Seevers, 1966, Kenyon, 1997, Vilhelmsen et al., 2014). The factor c is site-specific; it is in the range of 0.003–0.03 m/s3 for a = 1 and in the range of 3 to 11 m/s3 for a = 4 (Mohnke and Yaramanci, 2008). Günther and Müller-Petke (2012) found a value of about 0.0055 m/s3 (a = 1) for the East Frisian Island of Borkum, which is geologically very similar to Langeoog. If comprehensive calibration data for estimating KMRS are not available, Eq. (1) can provide a reliable index for lithological characterization in saturated rock, as relative changes in KMRS can nevertheless be identified from analysing MRS data. (Siemon et al., 2015). Due to the natural heterogeneity of the sediments and the methodological limits of conductivity determination methods (e.g. the inherent simplifications and assumptions), values for hydraulic conductivity have to be considered to be approximations with a certain range of uncertainty (e.g. Tietje and Hennings, 1996). Consequently, significant changes in hydraulic conductivity are usually observed on a logarithmic scale and thus we define an MRS lithology index derived from Eq. (1) using its logarithmic form:
According to Kenyon (1997), the parameters a and b in this study are set to 4 and 2, respectively. In principle, the parameter c can arbitrarily be chosen. It is convenient to determine c such that LMRS is always positive. In our case we chose 104.5. It must be noted that LMRS can only be used to characterize the saturated zone, because the decreased water content in the unsaturated zone must lead to an underestimation of LMRS. A second limitation of using LMRS for lithological characterization is that LMRS is not able to distinguish between lithological units, if these exhibit T2⁎ times smaller than the dead time of the MRS measurements. They will appear with similar signature.
For most of the MRS measurements on the island of Langeoog in 2012 and 2013, square loops with 50 m side length (2 turns) were used, mostly in combination with one or two remote reference loops for noise reduction (Walsh, 2008, Müller-Petke and Costabel, 2014). At three sites the noise was too high using this field setup, even after processing using the remote references. At these sites, figure-of-eight loops (with diameters 30, 40 and 50 m depending on the available space) for both measurement and reference loops were chosen (Trushkin et al., 1994, Müller-Petke and Costabel, 2014). In this way, an acceptable noise level could be achieved, but at the expense of decreased penetration depth compared to the normal setup (Trushkin et al., 1994). The signal length was set to 1 s for all readings. The post-processing for the MRS data consisted of several steps according to the procedure proposed by Costabel and Müller-Petke (2014): despiking, harmonic noise compensation, stacking, and envelope identification. Afterwards, the data were resampled using 50 logarithmically-distributed gates. This number of gates is larger than the sufficient minimum suggested by Dalgaard et al. (2016). However, a larger number of gates may stabilize the inversion calculation, because possible unwanted filter effects and artefacts are attenuated. For each gate an individual error is calculated and taken into account in the inversion (Günther and Müller-Petke, 2012). We used the Software MRSMatlab (Müller-Petke et al., 2016) for processing and inversion. The sensitivity functions were calculated on the basis of the resistivity distributions resulting from the TEM measurements conducted at the MRS sites. The instrumental dead time with the used equipment was 40 ms, while the effective dead time varied between 60 and 80 s, depending on the specific measurement and post-processing parameters that were adjusted individually for each measurement position to account for their individual noise characteristics.
For inverting the MRS data, we applied a smoothness-constrained 1-D inversion assuming a mono-exponential behaviour of each layer. The degree of smoothness of the resulting models is controlled by the regularization factor λ (Müller-Petke et al., 2016), which was varied from 1000 to 1000,000 for every data set. Finally, an optimal λ was chosen manually for each data set depending on its individual noise level. Similar to the procedure of inverting the TEM data, the smooth inversion results were analysed to determine the number of layers (three or four layers) and to generate the starting models for a subsequent block inversion.
Relationship of electrical resistivity and water content
The electric resistivity ρbulk of water-saturated rocks is related to the porosity Φ and the electric resistivity of the pore fluid ρfl. Archie (1942) suggested the power law
with an empirical exponent m to describe the relationship between ρbulk, ρfl, and Φ. This equation is valid if the pore water is the dominating component contributing to the current flow, as for sand and sandstone. It is not valid if the surface conductivity of particles is significant, e.g. for clay-bearing sediments. However, in this study we analyse solely the ρfl-values of the sand aquifers, which exhibit negligible clay contents.
Eq. (3) demonstrates that it is not possible to define a general value for ρbulk representing a sharp threshold to separate fresh from saltwater-bearing sediments, because the formation factor is a material-dependent parameter. For this work, we chose a value of 30 Ωm. With a formation factor of F = 4 typical for loose sand, this corresponds to a groundwater conductivity of 1333 μS/cm, which is reliable when considering that typical upper limits for drinking water are found in the range around 2000 μS/cm (at 25 °C).
Results and discussion of the individual method
Overview based on HEM data
After inversion, the HEM provides an overview of the 3-D spatial distribution of the electrical resistivity ρ for the entire survey area. Fig. 3 shows the resistivity distribution at five different depths (from 0 to − 60 m asl). The freshwater, indicated by ρ > 30 Ωm (blue and green colours), appears as several independent lenses. The main bodies are located in the west and in the east of the island, while in-between a small lens occurs (Fig. 3a to c). Probably due to the continuous extraction of groundwater, the western lens is shallower than the (unused) lens in the east. In the west, only small spots with ρ > 30 Ωm can be found at − 20 m asl (Fig. 3c), while almost the whole lens in the east shows freshwater at this depth level. Areas coloured yellow and orange are difficult to characterise, as the corresponding resistivity range (3 to 20 Ωm) can be interpreted as brackish water in a sandy environment or clayey or silty layers containing freshwater. Consequently, we cannot identify the nature of the lower boundary of the western lens from Fig. 3c alone. It is either delimited by an aquitard or by a zone of brackish water. At − 45 m asl, in the western part of the island, freshwater is not present anymore except of a tiny spot in the north-west. Apart from that, we only note brackish and salt water (Fig. 3d). In contrast, a certain amount of freshwater can be found in the eastern lens at the same depth. At − 60 m asl (Fig. 3e), the pore space under the island is almost completely filled with saltwater (ρ < 1.5 Ωm).
1-D example data sets from TEM and HEM
We have chosen the measurement position 1 (Fig. 1b; Table 1) at a distance of about 100 m in the north of the borehole at position 2 (Fig. 1b; Fig. 2a) to demonstrate and compare the different methods in detail. Fig. 4 shows and compares inversion results of the corresponding TEM measurement with the mean of the nearest 11 HEM data points, which are located about 70 m west of position 1. We calculated different models by smooth (SI) and block inversions (BI) for both TEM and HEM to demonstrate the level of ambiguity. Two BI results for each method are presented, each with different starting model. One starting model was defined according to the working flow as described above, the other one was defined with an additional layer representing an aquitard structure with decreased resistivity. In this way we want to analyse the impact of the starting model on the final result.
All resulting resistivity distributions in Fig. 4 start with a resistivity of about 100 Ωm at the top (sand with freshwater) and decrease to 1 Ωm at a depth of about 50 m (sand with saltwater). Near the surface, the HEM results show a small layer with increased resistivity (more than 1000 Ωm) and thicknesses of 2 to 3 m, which was included as pseudo-layer in the inversion to avoid possible systematic errors by disturbed measurements of flight altitude (Siemon et al., 2009). The SI of the TEM data shows an increase of the resistivity to 200 Ωm at the depth range of 10 to 20 m. We interpret this to be an artefact, because the resolution of the TEM method at shallow depths below 20 m is too low for resistive structures of > 100 Ωm. Except for the two features described above, the SI results of TEM and HEM are in general agreement. Both show a slight decrease of the resistivities down to 10–20 Ωm at a depth of approximately 25 m and a more or less homogeneous resistivity down to 40 m. Below 40 m, the resistivity decreases monotonously to about 1 Ωm at 50 m depth, which represents the transition zone to the saltwater.
A similar trend is shown by the BI results. For these, the transition zone between the freshwater and the saltwater bearing sand is represented by one or two single layers with a resistivity of 10 to 30 Ωm. The thicknesses of these layers vary depending on the individual starting model of each inversion. According to the lithology information of the borehole at position 2 nearby (Fig. 2a) and supported by the MRS data measured at the same position (discussed below), a clay layer at about 20 m depth with a thickness of several meters is assumed at position 1. However, neither the TEM nor the HEM measurement can resolve this clay layer unambiguously. On the other hand, as demonstrated by the five-layer TEM model, the assumed presence of a clay layer is not contradictory to the measurements. This model fits the measured data with the same RMS value as without the clay layer. However, this layer with decreased resistivity occurs in the inversion result only if the starting model already suggests a corresponding layer at the correct depth range.
MRS data with clay indications
The MRS inversion models (SI and BI) at position 1 are shown in Fig. 5, consisting of the corresponding distributions of θMRS (Fig. 5a) and T2⁎ (Fig. 5b) with depth. The SI result shows a θMRS of approximately 20% at shallow depth interpreted as the unsaturated zone. The corresponding T2⁎ is 0.3 s. In the saturated zone at z > 3 m, θMRS increases to about 28–35%, while T2⁎ decreases to about 0.2–0.25 s, a range typical for sandy aquifers. At about 30 m depth, θMRS decreases to 20%, whereas the T2⁎ distribution does not show a significant change. We interpret this with the presence of a clay layer, causing the MRS signal to relax within the measurement dead time (e.g. Boucher et al., 2011, Dlugosch et al., 2011). Compared to the lithological profile of the borehole at position 2 (Fig. 2a), where the clay layer occurs at a depth of 22 to 28 m, the clay indication by MRS is found deeper. In the second aquifer below the clay, which mainly consists of coarse sand, T2⁎ increases to about 0.5 s, while θMRS increases to 38%.
In accordance to the general trend of the SI result, the number of layers for the BI model was set to four: unsaturated zone, upper aquifer, clay layer, lower aquifer. The final result of the BI is in accordance with the SI result when taking the model uncertainties into account (red lines). The large BI model uncertainties for the water content of the first and the third layer demonstrate that the MRS data does not contain enough reliable information in the corresponding depth ranges.
Figs. 5c and d show the distributions of LMRS as derived from the SI and the BI result, respectively. The upper aquifer consisting of fine sand at a depth of 3 to 26 m is characterised by LMRS values between 1 and 1.5. The clay is indicated at a depth of 26 to 33 m by LMRS values smaller than 1. For the lower aquifer beneath the clay layer, mainly consisting of coarse sand and gravel, LMRS is found to be in the range of 2 to 2.5. SI and BI distributions of LMRS are in general agreement. The SI result exhibits a thin region with decreased LMRS at about 10 m. However, the lithology profile of the nearest borehole (position 2) does not show any lithology changes at shallow depths (Fig. 2a), so the origin of this feature remains unknown. Table 2 summarizes the relations between LMRS and lithology to present a key for the MRS data interpretation in the following analyses. In the unsaturated zone (down to 3 m in Fig. 5c and d) the decreased water content causes an underestimation of LMRS with values similar to clay as expected beforehand, which is omitted in Table 2 to avoid misinterpretations.
MRS data without clay indication
In contrast to the MRS data with a clay indication as discussed above, we now show an example without any aquitard indications. Fig. 6 shows the SI and BI results of an MRS measurement at position 3 (Fig. 1b; Table 1). As for the example at position 1, both SI and BI are in principle agreement considering the BI uncertainties. At shallow depths between 0 and 15 m, the θMRS and T2⁎ distributions start at about 30% and 0.2 s, respectively. At depths between 15 to 30 m, θMRS increases, while T2⁎ does not change significantly. The corresponding LMRS varies between 1 and 1.5 and indicates the fine sand aquifer as for the data example at position 1. At depths larger than 30 m, T2⁎ increases up to 0.7 s with θMRS decreasing down to 25% resulting in an LMRS of 2 for the SI result, whereas the LMRS of the BI remains at values below 2. The nearest borehole to this MRS measurement is found at position 4 (Fig. 1b; Table 1) at a distance of about 50 m in the south, which shows fine sand down to a depth of 20 m followed by medium sand down to 30 m. Coarse sand is found at a depth of 30 m down to the final depth of 52 m. We attempted to enforce the existence of a thin clay layer in the BI result by including a corresponding layer with typical NMR parameters in the starting model. However, the assumed clay layer did not appear in the final result. Instead, the inversion solution converged to a model where the included layer exhibited similar parameters as layer 1 and 2 in Fig. 6. Obviously, the presence of an aquitard layer cannot be enforced by a prejudged starting model, if not indicated by the MRS data.
Accuracy of clay indication depending on noise level
The data quality of the MRS data strongly varies in the area of investigation due to the variable noise conditions. Under optimum circumstances, the noise level was approximately 50 nV, but ranged up to 1000 nV near the power-lines connecting the production wells. To evaluate how the noise conditions affect the potential of the MRS method to indicate thin aquitard layers, a numerical study based on synthetic models was performed. Two synthetic block models and the corresponding results of simulated MRS measurements are depicted in Figs. 7a and b. The first model (Fig. 7a) consists of three layers, which represent a simplified geology: a coarse sand aquifer with θMRS = 0.4 and T2⁎ = 0.5 s overlain by a thin clay layer of 6 m thickness with θMRS = 0.45 and T2⁎ = 0.005 s, while the upper find sand aquifer is represented by θMRS = 0.35 and T2⁎ = 0.2 s. The second model represents the two-layer situation where the aquitard is absent. Using the same sensitivity function and dead time of the MRS measurement at position 1, a forward modelling was performed with 24 pulse moments in the range of 0.1 to 5.1 As. This data set was superposed by Gaussian-distributed noise of varying amplitude (10, 100, 500, and 1000 nV) and inverted afterwards. The resulting LMRS-distributions for the first model shown in Fig. 7a demonstrate that even for the highest noise level an indication for the clay layer can be found. The boundaries of the clay layer become increasingly blurred at higher noise levels, a trend that was expected beforehand. Additionally, with increasing noise level, i.e. increasing smoothness, the corresponding depth range with LMRS < 1 is overestimated more and more compared to the original model. The simulations for the second model (Fig. 7b) demonstrate that high noise levels do not generate artefacts, which might be misinterpreted as clay indications. Instead, the smearing of the boundary between both layers increases with the noise level as expected, while the LMRS in the transition zone never drops to values lower than 1.
Four MRS and four TEM measurements were arranged on a cross section in the vicinity of the HEM flight line FL16 (see Fig. 1b and Table 1), heading northwards. Both the SI (Fig. 8a) and BI (Fig. 8b) results of HEM and TEM are in good agreement. The freshwater lens (ρ > 30 Ωm, green to light blue colours) is clearly visible from 100 to 700 m on the profile and reaches a mean elevation of − 10 to − 20 m asl. The saltwater (ρ ~ 1 Ωm, red to purple colours) is found at − 40 to − 50 m asl. Between freshwater and saltwater, we note a transition zone with resistivities ranging from about 30 down to 1 Ωm for the SI. The thickness of the transition zone could be verified at the borehole at position 18, which is discussed below. For the BI, a layer with resistivities ranging from 10 to 20 Ωm (yellow colours) indicates the transition zone. In the northern part of the profile, both the SI and the BI profiles show a layer with a decreasing resistivity (5 Ωm and less) at − 10 to − 20 m asl. Towards the Wadden Sea in the North, this layer rises monotonously, while its resistivity decreases down to the level of saltwater-saturated sediments. Due to the ambiguity when interpreting resistivity information, the nature of this layer cannot be identified by the EM methods alone. Additional information is needed and provided by MRS.
Fig. 8c shows the LMRS results of the MRS measurements at the positions 3, 5, 6, and 9 together with the lithology of the boreholes near FL16 (from south to north: 10, 8, 4, and 2). Based on LMRS and the borehole information, two separate aquitards can be identified, one in the north and one in the south of the profile. Their specific lithology is known from the drilling logs of the boreholes, i.e. clay in the North (position 2, see Fig. 2a) and a sequence of fine sand, silt, and clay in the South (position 10, see Fig. 2b). The spatial dimensions of the aquitards can be reconstructed with the additional information from LMRS. Moreover, the LMRS distributions allow distinguishing between the fine sand at shallow depths down to about − 20 m asl and medium to coarse sand at − 30 to − 40 m and below.
When plotting the reconstructed aquitard layers onto the HEM section in Fig. 8a (black dashes), we note that in the north the layer with decreased resistivity is on top of the clay layer, while the clay layer itself is not clearly outlined. Our interpretation is that brackish or salty water being trapped on top of the clay, a situation that is frequently seen at coastal aquifers (Günther and Müller-Petke, 2012, Attwa et al., 2011). As the measured resistivities above the clay layer decrease with increasing distance from the sea shore, we believe that the brackish water on the clay is diluted and/or displaced more and more by infiltrating freshwater. Another hydro-dynamically relevant effect is observed in the middle of the profile: in the gap between the northern and the southern aquitard, i.e. between 250 and 400 m on the profile line, the thickness of the transition zone between freshwater and saltwater is smaller compared to the transition zone beneath the aquitard. Without the aquitard, infiltrating freshwater can move faster downwards and thus, the freshwater lens reaches a greater depth (− 30 m asl).
On the profile along the flight line FL17, five TEM measurements, three MRS measurements, and two boreholes are available (Fig. 9). Both the SI and the BI results of the HEM data show that the lateral dimension of the freshwater lens (170 to 720 m on the profile) is smaller compared to FL16, while its mean thickness (30 to 35 m) is significantly higher. The transition zone to the saltwater seems to be absent or at least very thin in the middle of the profile between 200 and 400 m. Outside this range, a transition zone of similar thickness as on FL16 is present. Again, the results of HEM and TEM are in good agreement.
In the middle of the profile between 250 and 500 m, the MRS measurements and the borehole at position 12 indicate an isolated aquitard that is not indicated by the EM measurements. On the other hand, an aquitard in the north of the profile is indicated only by the HEM results: a layer of decreased resistivity similar to the northern structure in FL16. In accordance to the interpretation of FL16, we assume that also on FL17 salt/brackish water is trapped here on top of a clay layer. TEM/MRS measurements were not possible in this area due to its steep and undulated topography, thus no estimates of the thickness for the clay layer are available.
The only E-W flight line in the area of investigation is FL2.9 (Fig. 1b). Its SI result is depicted in Fig. 10 together with the lithology of all boreholes at a distance of less than 100 m from the line. In addition to the borehole information, the results of three MRS measurements in the vicinity of FL2.9 are included (positions 3, 14, and 15 in Fig. 1b). Both MRS results, i.e. the corresponding distributions of LMRS, are already depicted in Figs. 8c and 9c. Thus, Fig. 10 shows the lithological interpretation of the MRS measurements using the same format as for the borehole lithology logs. As for FL16, we also note that freshwater can be found down to − 20 m asl (from 500 to 3200 m on the profile line). At two positions (x ~ 1600 and 2200 m) and in a larger region in the eastern part (x = 2400 to 3000 m), the freshwater reaches elevations of about − 30 to − 35 m asl. This increase in thickness was also observed on FL17 (Fig. 9). A layer with decreased resistivity (< 5 Ωm) is observed at about − 20 to − 30 m asl from x = 1000 to 1500 m and from x = 1700 to 2100 m, which indicates a clayey aquitard, verified by the lithology logs of observation boreholes 27 and 25. The HEM profile clearly indicates the interruptions of the aquitard, e.g. at x = 1600 and 2200 m. In the east, the HEM data suggest that the aquitard is absent. Although this observation is verified by the borehole and the MRS measurement at positions 4 and 3, the borehole 2 and the two eastern MRS measurements show an aquitard where the HEM does not. However, the borehole at position 2 is located 75 m away from the HEM line, so the lithology is not necessarily the same. Regarding the MRS-based aquitard indication at 14 and 15, we note again as in Fig. 8, Fig. 9 that the aquitard structure consisting of alternating sequences of fine sand, silt, and clay is not visible in the HEM image.
Mapping of aquitard layers
Using the lithological information from boreholes and MRS measurements, a map showing the aquitard indications is given in Fig. 11. As described above, the MRS measurements are concentrated in the northeast of the investigation area, where the borehole density is very sparse. The MRS-based information about the aquitard occurrences complements the existing geological database of the managed freshwater lens and can be implemented in the hydrogeological model. In addition, the HEM data of the E-W flight line FL2.9 (see Fig. 10) strongly suggests the existence of an aquitard between about − 25 and − 30 m asl also in the centre of the area. This was unknown to date, because the depths of the boreholes in this area are less than 20 m. However, the new geophysical insights from HEM and MRS verify the findings and expectations of previous studies (e.g. Tronicke, 1997, Houben et al., 2014). Although the extent of aquitard layers seems to be larger than expected beforehand, large gaps allow rapid water flow from the shallower fine sand aquifer to the deeper medium sand. Although of less relevance for groundwater modelling at regional scale, the knowledge of these gaps for local water suppliers is of enormous importance when assessing for instance the risk of salt water upconing at specific production wells.
Estimation of the electrical conductivity of groundwater from geophysical measurements
De Louw et al. (2011) and Gunnink et al. (2012) showed that the SI results of HEM are in good agreement with electrical cone penetration tests. Thus, we assume that our EM results provide reliable estimates of the formation resistivity and can be used to predict the electrical conductivity σfl of the groundwater, i.e. the reciprocal of ρfl (Eq. (3)). Using the ρ results of HEM and TEM together with θMRS as rough estimates of Φ, σfl can be estimated from the geophysical data using Archie’s equation (Eq. (3)). To do this, we set the empirical parameter m to 1.3, a reference value for sand that was also determined at the neighbouring island Borkum by using direct push and fluid conductivity measurements (Günther and Müller-Petke, 2012). Figs. 12a and b show the SI results of TEM, HEM, and MRS for the measurement point 6. This point is located 70 m southeast of a borehole (position 18, see Fig. 1b and Table 1) with an unusually long screen covering a depth between 35 and 65 m, thus providing the opportunity of observing the transition zone between fresh and saltwater directly. The corresponding MRS (Fig. 12 b) indicates the existence of a clay layer as found at position 1 in Fig. 5. The clay is indicated by decreased water content at depths shallower than the screen of borehole 18, i.e. at the depth range of the screen the requirements of Archie’s law (Eq. (3)) for estimating σfl are satisfied. Thus, this borehole provides ground truth to demonstrate the estimation of the vertical σfl distribution using the EM and MRS data. The fluid conductivity in borehole 18 was measured using a salinity logging tool. As it was most likely uninspected for decades, the water salinity within the borehole was expected to approximate the σfl of the formation. After the logging, three depth-specific groundwater samples were taken for verification, each after 1 hour of pumping. Assuming a primarily horizontal flow, the corresponding σfl of these samples are expected to represent the original groundwater values. Estimates and measurements of σfl within the depth range of the screen depicted in Fig. 12c are in good agreement. At the upper parts of and above the screen, the σfl graph representing the logging data shows smaller σfl values compared to the reference measurements after sampling and compared to the geophysical estimates. We assume that this is a result of diffusion and dispersion of dissolved ions inside the borehole.
In addition to the observation at borehole 18, further groundwater samples were investigated by sampling of observation boreholes (Houben et al., 2014). All these samples are taken from short screens in the sand aquifers at different depths covering a broad range of σfl in the area of investigation. We use the σfl data at boreholes next to the TEM/MRS measurement positions to test how accurate σfl can be estimated using the geophysical data. Fig. 13a correlates the measured σfl with the σbulk as measured by HEM and TEM. The obvious linear shift between both quantities suggests that using the traditional approach of considering a constant formation factor would lead to an accurate estimation of σfl from the EM data, once the formation factor F (Eq. (3)) is determined. Estimates of F for sediments and rocks can be found in the literature. However, a correct determination requires a comprehensive study of the σfl-σbulk-relationship in the whole area of investigation. Thus, we also test how accurate the fluid conductivity can be estimated using the resistivity data from the EM methods and the water content from the MRS data. In this case only the Archie-exponent m (Eq. (3)) must be determined, the variability range of which is expected to be much smaller than that of F for loose sediments. Fig. 13b shows the corresponding cross-plot for the estimates from the combination of HEM/MRS and TEM/MRS. All estimates match the actual values within half a decade and with a coefficient of correlation (R2) of 0.98 for TEM and 0.95 for HEM, which suggests that including the MRS data for the estimation yields reliable σfl estimates without the need of calibrating the formation factor.
Fig. 11 and the W-E profile in Fig. 10 show that aquitard structures with thicknesses of up to 10 m at about − 15 to − 20 m asl seem to be common in the entire area of investigation. In the east of the area, where boreholes are sparse, the additional information provided by MRS complements the few available lithological logs. However, the thicknesses of the aquitard layers estimated by MRS must be handled with care. As the numerical study demonstrates (Fig. 7), aquitard indications from noisy MRS data exhibit overestimated thicknesses, even for low noise levels. Aquitard indications were also found in the EM data. As demonstrated in Fig. 10 (HEM) and in Fig. 8, Fig. 9 (TEM and HEM), clayey layers can show up through decreased resistivity either by their lithology or by salt/brackish water trapped on them. In the second case it is not possible to identify the aquitard itself by EM. In such cases the combination of EM with MRS is very promising. However, distinguishing between clay and silt aquitards was not possible with MRS in this study, because both exhibit relaxation times smaller than the effective dead time. Recent technical developments have led to much shorter instrumental dead times down to 4 ms (Walsh et al., 2011). However, the effective dead time is longer due to relaxation during pulse and increases with the pulse length (Walbrecker et al., 2009, Dlugosch et al., 2011). Thus, the problem of characterizing water in small pores remains an issue for deep groundwater exploration at depths of several deca-metres, because long pulses are necessary to reach these depths. For shallower depths below 30 m, MRS measurements with instrumental dead times shorter than those applied in this work can identify groundwater also in silt formations as demonstrated by Günther and Müller-Petke (2012).
In spite of the advantages provided by MRS regarding its sensitivity to lithological changes in the subsurface, a serious shortcoming must be noted. Due to the high noise sensitivity of the MRS method and the corresponding high stacking rates required, the measurements are much more time-consuming compared to EM methods. Only one to two MRS measurements per day could be performed. It is not expected that the technical equipment can be improved significantly in the near future to overcome this problem. Thus, future hydrogeological investigation programs on the regional scale with MRS involved cannot be conducted using large MRS profiles or even 3-D observations with an adequate depth resolution. Instead, the usage of MRS will be limited to single point measurements with high vertical resolution or, on the other hand, to profile measurements with reduced vertical resolution (Costabel et al., 2016). The latter is only feasible if the subsurface can be described as a layered model. The preferred solution in general could be the combined use of MRS with HEM: the regional overview is given by the 3-D resistivity distribution from HEM, while MRS (in addition to borehole information) helps to overcome the natural ambiguity of resistivity data regarding the lithological interpretation at the local scale, i.e. at specific places of importance.
The MRS-based lithology index LMRS that we derived from the well-known empirical relationship between hydraulic conductivity and NMR data (Eqs. (1), (2)) does a good job in combining both NMR parameters θMRS and T2⁎ in one handy measure and to link them to lithology. As for the usual empirical quantification of hydraulic conductivity from NMR and MRS data, we must expect the relation of LMRS and lithology to be site-specific. However, a recent study of Knight et al. (2016) showed that the empirical parameters of Eq. (1) differ only little in unconsolidated sediments for 10 different investigation sites. Consequently, the assumption is reasonable that LMRS can be treated as a general measure for lithology without site-specific “calibration”. Therefore, future research should investigate the variation of LMRS for different sites and geological settings.
In this case study, we applied the geophysical methods of HEM, TEM, and MRS, complemented by lithological and hydrochemical data, to investigate the freshwater lens on the North Sea island of Langeoog. Besides the dimensions of the lens, we extracted lithological information from the geophysical data to complement the sparsely available borehole information. Of specific interest was the spatial distribution of aquitard layers that have a great impact on the internal dynamics of the freshwater lens. We found that especially MRS is able to provide reliable indications for the occurrence of clay and silty fine sand layers with thicknesses of less than 10 m.
We conclude that a combined application of EM and MRS helps resolving the ambiguity between lithology and salinity effects from the resistivity data. Additionally, reliable estimates of groundwater electrical conductivity are provided by the EM and MRS data for the sandy aquifers. Future improvement of the data interpretation might be achieved by jointly inverting EM and MRS data (Behroozmand et al., 2012, Günther and Müller-Petke, 2012). As HEM is appropriate to close the gaps between the often sparse information of ground-based MRS data, strategies to combine the 1-D and 3-D data sets in one inversion scheme are necessary, similar to the approaches of Gunnink and Siemon (2014) who included lithological information from boreholes into the HEM inversion, or Sudha et al. (2014), who combined TEM and HEM in a joint inversion algorithm. Our analyses show that the results of HEM and TEM are very similar, which suggests that the HEM resistivities are suitable to be used for the MRS sensitivity calculations without the additional effort of conducting TEM measurements in the future.
The lithological data extracted from the geophysical data set demonstrated in this study will be included in a future numerical groundwater flow model for the island of Langeoog. Filling the gaps between often sparse geological data from boreholes using additional geophysical information will help improving the reliability of corresponding simulations. This will support the water supplier in managing the fragile freshwater resources properly and in planning specific actions for future scenarios.
Source: Geophysical investigation of a freshwater lens on the island of Langeoog, Germany – Insights from combined HEM, TEM and MRS data
Authors: Stephan Costabel, Bernhard Siemon, Georg Houben, Thomas Günther