Hydraulic properties at the North Sea island of Borkum derived from joint inversion of magnetic resonance and electrical resistivity soundings

For reliably predicting the impact of climate changes on salt/freshwater systems below barrier islands, a long-term hydraulic modelling is inevitable. As input we need the parameters porosity, salinity and hydraulic conductivity at the catchment scale, preferably non-invasively acquired with geophysical methods. We present a methodology to retrieve the searched parameters and a lithological interpretation by the joint analysis of magnetic resonance soundings (MRS) and vertical electric soundings (VES). Both data sets are jointly inverted for resistivity, water content and decay time using a joint inversion scheme. Coupling is accomplished by common layer thicknesses. We show the results of three soundings measured on the eastern part of the North Sea island of Borkum. Pumping test data is used to calibrate the petrophysical relationship for the local conditions in order to estimate permeability from nuclear magnetic resonance (NMR) data. Salinity is retrieved from water content and resistivity using a modified Archie equation calibrated by local samples. As a result we are able to predict porosity, salinity and hydraulic conductivities of the aquifers, including their uncertainties. The joint inversion significantly improves the reliability of the results. Verification is given by comparison with a borehole. A sounding in the flooding area demonstrates that only the combined inversion provides a correct subsurface model. Thanks to the joint application, we are able to distinguish fluid conductivity from lithology and provide reliable hydraulic parameters as shown by uncertainty analysis. These findings can finally be used to build groundwater flow models for simulating climate changes. This includes the improved geometry and lithological attribution, and also the parameters and their uncertainties.


Introduction
Climate changes are threatening the groundwater dynamics, particularly at the coast and on barrier islands, where the interaction between freshwater and saltwater is important for water supply.This is the subject of the international and interdisciplinary project CLIWAT (CLImate and WATer) investigating the impact of sea-level rise on freshwater resources at the coasts of the North Sea and Baltic Sea.In one of the projects, modelling of the long-term hydraulic behaviour of the freshwater lens beneath the island of Borkum is accomplished (Sulzbacher et al., 2012).Besides lithological models, for density driven flow models several input parameters are required such as porosity, hydraulic conductivity and salinity.Very often these quantities are not available at the catchment scale.
Point information can be retrieved by boreholes and other direct investigations.For large modelling, however, these are usually not available sufficiently dense.Geophysics can play an important role in closing the gaps between boreholes.Airborne electromagnetic measurements are particularly important since they provide three-dimensional distribution of electrical conductivity, or its inverse, resistivity (Siemon et al., 2009), which is a key parameter in hydrogeophysical investigations (e.g.Viezzoli et al., 2010).However, there are two main shortcomings of resistivity: (i) we cannot clearly differentiate between clay content and fluid salinity, and (ii) there is no sufficiently reliable relation to hydraulic conductivity, probably the most important parameter needed.
Geophysical techniques based on the principles of Nuclear Magnetic Resonance (NMR) can contribute to overcome these shortcomings.The method measures a signal released by set of hydrogen protons that relax from an excited state back to equilibrium.This relaxation process can be described by exponentially decaying functions.NMR allows for uniquely determining water content of a sample based on the direct sensitivity on the number of hydrogen protons represented by the initial amplitude of the exponential function.Furthermore, the measured decay time depends on the pore geometry and can therefore be used to estimate permeabilities (Seevers, 1966).There are several types of relaxation and corresponding experiments.Whereas both longitudinal relaxation time T 1 and transversal relaxation time T 2 mainly reflect pore geometry, the free induction decay (FID) time T * 2 is additionally affected by magnetic gradients.
NMR is well known and established as a laboratory and borehole method and with increasing success applied at the field scale.Surface NMR experiments utilize large surface loops to produce electromagnetic pulses with increasing pulse moments q (product of current and duration), which successively reach deeper parts of the subsurface.For a detailed explanation of the method see, for example, Legchenko and Valla (2002) and Hertrich (2008).Particularly for parametrisation of hydrogeological systems see Lubczynski and Roy (2004) and Lachassagne et al. (2005) and references therein.
A coincident loop experiment is referred to as magnetic resonance sounding (MRS), since inversion retrieves water content and decay time in the subsurface as a function of depth.The most general approach was presented by Mueller-Petke and Yaramanci (2010) using the full data cube along the pulse moment (q) and time (t) axis.This class of inversion is therefore referred to as QT inversion.Mueller-Petke and Yaramanci (2010) discretise the subsurface in the spatial (depth) and spectral (decay time) dimension and achieve a smooth distribution of T * 2 for each of the depth layers.However, very often the subsurface consists of distinct layers of constant properties, and a mono-exponential decay is a valid assumption for many unconsolidated sediments (Hertrich, 2008).Therefore, we follow a block scheme inverting for the parameters water content and decay time of a small number (typically between two and six) of layers with variable thickness.This approach strictly uses monoexponential decay compared to the scheme with stretchedexponential decay very recently presented by Behroozmand et al. (2012a).We argue for using the simplest model that satisfies the data, in our case mono-exponential behaviour showed to be sufficient.
For calculation of MRS responses, a resistivity model is needed to determine the magnetic fields in the subsurface (Weichman et al., 2000).Theoretically, the resistivity could be retrieved from the phase information of the measurements themselves (Braun and Yaramanci, 2008).A simultaneous inversion for the three parameters, i.e. water content, decay time and resistivity, was done by Braun et al. (2009) using a time step inversion approach and a distinct number of layers with constant water content and mono-exponential decay times.However, only resistivities below 25 m effectively influence the signal (Braun and Yaramanci, 2008) for a loop diameter of 50 m.Depending on the geology this can be sufficient to distinguish aquifer and aquiclude structures, but it is unlikely to distinguish unsaturated from saturated sand.Therefore, we prefer to combine MRS with a direct current (DC), frequency domain electromagnetic (FDEM) or transient electromagnetic (TEM) sounding.For resistivity interpretation a block discretization is typical and can be solved by linear filtering and fast Hankel transformations (e.g.Anderson, 1989).
Since both methods are sensitive to the typical structures of an aquifer system, a combined or joint inversion is favourable.The coupling of the methods is achieved only by the common layer thickness.Hertrich and Yaramanci (2002) presented a joint inversion scheme for resistivity and water content using a generalized Archie model and differentiated between bound and mobile water.Our objective is to further include the decay times of NMR signals for both structural identification and hydrological characterisation.
The relations to obtain hydraulic conductivity K from laboratory NMR measurements go back to the model of Seevers (1966) using a second order dependence of K on the decay times T 2 or T 1 and a first order dependence on the porosity .There are only a few papers on retrieving hydraulic conductivity K from free induction decay (T * 2 ) measurements in the field scale.For an overview see, for example, Mohnke and Yaramanci (2008) and Plata and Rubio (2008) and references therein.T * 2 is sensitive to inhomogeneities of the magnetic field, i.e. magnetic gradients at the pore scale reduce T * 2 (Grunewald and Knight, 2011).Therefore, most MRS papers deal with T 1 relaxation times to avoid this problem (Legchenko et al., 2004).However, Walbrecker et al. (2011a) demonstrated that even a small off-resonance excitation (which is inevitable in field experiments) has large effects on the double-pulse experiments used to retrieve T 1 so far.These off-resonance effects can lead to significantly wrong parameters.Moreover, Walbrecker et al. (2011b) showed that the inversion scheme for T 1 measurements is not generally valid.Since an appropriate measuring scheme was not available at the time of our study, we did not use T 1 , but T * 2 measurements.Concerning the interpretation we need to be aware of the ambiguity of T * 2 , which is controlled by both magnetic field inhomogeneities and pore size.
Close to our work, Vouillamoz et al. (2007) jointly interpreted MRS and vertical electric soundings (VES) and characterized aquifers based on NMR parameters and resistivity.The authors showed uncertainties of the derived parameters and demonstrated that inverting VES with fixed geometry from MRS significantly improves resistivity uncertainty.Very recently, Behroozmand et al. (2012b) presented a study on joint inversion of MRS and TEM that is similar to our approach.While describing geophysical aspects and inversion schemes in detail, they do not discuss the derived quantities in the context of hydrological parameters.
In our paper we will present a methodology to invert for NMR parameters and resistivity simultaneously, applying a QT block inversion.After briefly presenting the methodology, we will show the results of three soundings measured on Borkum Island.One of them is used for verification of the method by comparison with a borehole.Another one is used for petrophysical calibration using a pumping test.Finally, we demonstrate on the third measurements how hydraulic quantities are retrieved and how big their uncertainties are.

MRS block QT modelling
We assume the resistivity model of the subsurface is known, e.g. from electric or resistivity soundings.The complex forward response (initial amplitudes) for a fixed discretisation can be formulated in terms of a matrix-vector multiplication: where ũ is the complex vector of simulated voltages, w is the searched water content vector and K is the complexvalued kernel.The latter depends on loop geometry and the resistivity distribution.Details about the computation can be obtained from Weichman et al. (2000) or Hertrich (2008).
Usually instead of the complex data, amplitudes are inverted for by transforming Eq. ( 1) into a real-valued matrix-vector equation, u = K w.The problem becomes non-linear, i.e.K depends on w.Since the kernel computation is extensive (calculation of magnetic field intensity B and integration from 3-D to 1-D), it is done only once for a relatively fine discretization.The forward response of an arbitrarily located layer is then derived by summing the original layers with their weight of coincidence between two models (Behroozmand et al., 2012a).QT type inversions (Mueller-Petke and Yaramanci, 2010) use the whole data cube along both the pulse moment (q) and time (t) axis.Each value of the response for a single layer, i.e. the amplitudes only from that layer, is then multiplied by the exponential function e −t/T , forming the data cube.This is done for all layers, and the values sum up yielding the forward response vector f , i.e. a data cube of amplitudes E for a number Q of pulse moments and a number G of sampling times: for a given model vector m.The latter contains, for a number of L layers, thickness values z i (i=1,. . .,L-1), water contents θ i (i=1,. . .,L), and decay times T i (i=1,. . .,L).
Even though the number of unknown parameters is very small (3 L-1), the number of data (G • Q) and thus the size of the inverse problem can become large due to the high sampling rate (typically 10 kHz) of the measured signals.Therefore, we follow Behroozmand et al. (2012a) and resample the individual decays into a number of about 40 gates using an integration procedure.Logarithmically equidistant gate lengths are defined and the arithmetic mean of all values within the gate is derived.Thus, after gate integration the statistical error of a gate value changes from gate to gate as a function of the individual gate length.We calculate this error by dividing the stacking error of raw data by the square root of the number of values being averaged within a gate.This is equivalent to the usual stacking improvement assuming Gaussian distributed noise.Each individual gate error is then taken into account (cf.Eq. 2) using error weighted data inversion.
For the sake of clarity, the data error before gating is referred to as the noise level and can be obtained from (i) pure noise measurements, (ii) from the residual imaginary part of the data after rotation, or (iii) from stacking (Müller-Petke et al., 2011).We decided for the latter since pure noise measurements were not available with the used instrument and it is the most general approach.

Inversion scheme
As typical for block inversion, we use a Marquardt-type damped Gauss-Newton inversion (Inman, 1975), i.e. using a local damping with successive cooling of the regularisation parameter.In order to account for the different measured quantities (apparent resistivities ρ a in m, and E(q, t) in V), we apply a data weighting using independent errors i for each of the total number of N data points (d i ) so that the objective function to be minimized reads The data covariance matrix C d contains the variances 2 i on the main diagonal.The data-weighted misfit can also be expressed as chi-square function, χ 2 = /N.A value of 1 means fitting the data within error bounds in a least-squares sense.In each inversion step, the update m to the model m is retrieved by solving where f denotes the forward operator, J is the Jacobian matrix, and I is the identity matrix.The regularization parameter λ is successively decreased until convergence is reached.Since the number of unknowns is very small, a derivative of the forward response (the sensitivity matrix) is easily obtained by the perturbation method, i.e. one additional forward run with a slight change for each individual value.
Regarding prior information on valid parameter ranges helps to improve the inversion results.Therefore, each of the 3 L-1 unknowns, p i , is transformed by a double-logarithmic transform (Kim et al., 1999;Günther, 2004) to the associated model parameter with p l i and p u i being lower and upper bounds of p i , respectively.The logarithmic transform has the advantage of automatically holding the values within bounds while decreasing the ill-posedness of the inverse problem.Table 1 gives an overview of the values used.They represent conservative bounds of values found in literature and could be further refined for each layer if prior knowledge is available.
The choice of starting models can be deciding since the algorithm can be trapped in local minima.In order to not drive the inversion too much, we use a homogeneous starting model of θ = 0.2 and T = 0.1 s.For resistivity the mean of the apparent resistivity is used.The value of the initial layer thickness turns out to be most critial.From the MRS kernel function we estimate a maximum depth comprising 80 % of the cumulative sensitivity (Christiansen and Auken, 2012) and divide this penetration depth by the number of layers.

Joint MRS and VES inversion
A joint inversion is straightforwardly achieved by concatenating all data quantities, i.e. the combined response vector becomes f = [f MRS f VES ] T .Hence, the Jacobian matrix obtains the correct form of two concatenated matrices.Accordingly, the variance vector for building C d is combined.Note that by the use of appropriate errors in the data covariance, a weighting of individual data fits is not needed any more.In all our inversions with χ 2 ≈ 1, we observed that the χ 2 values of both methods were close to 1 as well.
For a pure MRS inversion, we must assume a resistivity distribution.To avoid biasing the result by using a contrasted resistivity model, we decided to start with the multilayer inversion result of the closest helicopter electromagnetic (HEM) sounding (cf.Fig. 1), which was also used for the groundwater model of Sulzbacher et al. (2012).Since models are quite smooth, there is no initial support for layer boundaries that are equal to the resisitivity model.Whereas for the single MRS inversion this resistivity model remained constant, for a rigorous joint inversion we needed to couple the improved VES result with the MRS using a kernel update.As a consequence, the joint model changed, making an iterative approach necessary.Since a kernel update is expensive and to avoid the risk of being trapped in a local minimum, we updated the resistivity in an outer loop.Experiments showed that no more than 3 outer iterations, each fully minimizing the objective function, are needed.
Another important issue is the choice of the number of layers.In cases where no boreholes are available, we suggest to start with a homogeneous model and increase the number of layers successively until no further decrease in the objective function is observed.Finally, appropriate initial values might guide inversion but can also hinder convergence.

Computation of uncertainties
It is important for a parameter estimation to quantify the range in which the obtained values are expected to be.There are several ways of computing uncertainties: (i) based on the resolution matrix, (ii) based on the a-posteriori model covariance matrix, (iii) using most-squares inversion, and (iv) by a variation of individual parameters.As used by Müller-Petke et al. ( 2011), we decided to use the latter way, i.e. to vary the individual values independently unless the data confidence interval is exceeded.This method does not linearise the problem around the solution as other methods and yields different values for lower and upper bounds that are, however, close to the values derived from the model covariance matrix.After successful inversion each parameter is varied until the forward response deviates from the solution by an amount associated with χ 2 = 1.This corresponds in a statistical sense to a 95 % confidence interval.
Of course parameters are expected to infer each other.For example, the product of water content and thickness of a layer denotes the amount of water being proportional to the signal strength and is better described than both parameters alone.However, in a hydrological investigation this corresponds to a pumping test.Furthermore, the most critical parameter, the decay time, is expected to be relatively independent of the others.

Measurements and data processing
The eastern part of the divided freshwater lens of Borkum is the target of the investigation.In this area dunes of significant elevation prevail.These restrict the layout of large loops needed to reach the targeted investigation depth of about 60 m with MRS.Four soundings, of which three will be discussed, were acquired in spring 2010 within the frame of a BSc thesis (Liebau, 2010).The main objective of this survey was a feasibility study.Only amplitude inversions were conducted, and therefore only the locations of aquifers were derived but not their parameters.Figure 1 shows the location of the measuring and reference loops along with VES positions, HEM flight lines and the borehole.One of the soundings, CL2, was placed in the middle of the dune area, which describes the main extent of the freshwater lens.A research borehole (CLIWAT II) provides information about the lithology necessary for the validation of the method.Another sounding, OD33, was conducted at the southern boundary of the dunes.It is situated at a water well, where pumping tests have been carried out with the aim of calibrating the hydraulic conductivity equation.The last one, SKD, was placed at the eastern boundary of the freshwater lens where a significant silt layer was presumed.Table 2 summarises the acquisition parameters.
The geomagnetic field was about 49 300 nT and the corresponding Larmor frequency at about 2100 Hz.For all measurements, square loops with two turns (black rectangles in Fig. 1) have been used in order to increase signal strength.The loop dimensions were chosen to reach the target depths after sensitivity analysis, i.e. about 50 m edge length above the centre of the freshwater lens and about 25 m in the flooding area where the investigation depth was restricted by saltwater.We used the GMR instrument (Vista Clara Inc.), a multi-channel device with instrument dead times below 10 ms.The pulse lengths for the first three soundings were 40 ms in order to maximize the pulse moment to about 7 As, as necessary for a deep penetration.For the smaller target depth of the fourth sounding, we decreased the pulse length, which improves the ability of detecting fast decaying units such as silt (Dlugosch et al., 2011).
Additionally, reference loops have been laid out in order to register noise to be removed from the signals using noise cancellation.This technique makes use of a independent registration on another channel of a multi-channel instrument (Dlugosch et al., 2011).A transfer function is calculated to correlate the noise signal collected in the reference loop with the noise in the detection loop.This transferred noise signal from the receiver loop is then subtracted from the detection loop (Müller-Petke and Yaramanci, 2011).Doing so, a major part of the present noise at CL2 (about 80 %) was already removed prior to further processing.For OD33 and SDK no reference loops were needed since the noise level was already low.Generally, very low noise was observed with mean noise levels of 107 nV for CL2, 17 nV for OD33 (both after noise cancellation), and 9 nV for SKD (without noise cancellation).The noise level of measurements is an important input parameter in inversion (cf.Eq. 2), particularly if different data types are combined.
For the combined inversion the closest vertical electrical soundings have been chosen (see Fig. 1), having distances of not more than 50 m.The VES were all measured in the 1990s and again in 2008 (M.Grinat, personal communication, 2012).They comprise Schlumberger soundings with logarithmically increasing AB/2 starting from 1.5 m to about 150 m.Data quality was very good and the noise level was estimated with 3 % relative error.The latter is higher than the pure measurement error and includes electrode position errors and 3-D effects caused by small-scaled irregularities.

Verification sounding CL2
The first sounding (CL2) was conducted directly at the research borehole CLIWAT II, which was drilled in autumn 2009, six months before the MRS measurements.The borehole is located in the middle of the dune area where the freshwater lens reaches its maximum depth of more than Table 2. Acquisition parameters: loop sizes, number of pulse moments Q and stacks, maximum pulse q Q , pulse length τ p and effective dead time t e (see Dlugosch et al., 2011), maximum of the initial amplitudes max (E 0 ), noise after stacking, and S/N corresponding to max(E 0 ).

Name
Loop size 50 m.Lithology is well known from the interpretation of the drilling material.Furthermore, gamma ray and induction/resistivity log measurements were carried out directly after the drilling, having a hint about clay content and lithology (T.Wonik, personal communication, 2012).Figure 2 shows the different available resistivity models alongside the gamma log data and lithology interpretation.
On top there is a thick Holocene aquifer of clean, fine sand with the water table at about 3 m depth.It is followed by a silt-sand-clay layer representing the Holocene base.Below 32 m there is a second Pleistocene aquifer of brown-gray fine sand followed by an inter-bedding of sand and clay at about 50 m depth.The transition zone from fresh to saline water is also at this depth, as clearly mapped by the the array induction (AI) log or the vertical electrode chain (Grinat et al., 2010).The latter is a buried direct current resistivity instrument for continuous monitoring and represents the most reliable resistivity information since it is neither affected by borehole fluid nor by inversion ambiguity.All three surface soundings are generally able to detect the course of resistivity.However, they fail to yield a clear hint to lithology changes at this site.The ambiguity of possible models within data error is illustrated by the smooth 15-layer HEM model as used by Sulzbacher et al. (2012).
Resistivity variations in the freshwater regime (100 m or higher) do not significantly affect the MRS kernel (Braun and Yaramanci, 2008).From the borehole we have a very good resistivity model that made a kernel update obsolete.For joint inversion we chose a 5-layer model to account for the dry sand, the two aquifers, the aquitard and the conductive clay/saltwater zone.From the latter we do not expect to detect an NMR signal since decays from clay are too fast for current instruments and the high conductivity channels the magnetic fields thus reduce the amplitude of signals from below.Even though the data could probably be fitted equivalently by only four layers, we accounted for the available lithology information.
After 11 iterations the total error-weighted misfit decreased to χ 2 = 1.15.For both single and joint inversion, the data fit of the MRS and VES data reached similar values of about χ 2 = 1.Data fit and the results of the joint inversion are shown in Fig. 3.The uppermost layer is the unsaturated zone characterized by high resistivity and low water content.Otherwise, the water content variations are rather small (27-32 %).
The first aquifer is characterized by resistivities of about 100 m, 30 % water content and 200 ms decay time, which is relatively high for fine sand (Schirov et al., 1991;Lubczynski and Roy, 2003).Since the MRS measured T * 2 decay time is, besides the sensitivity to pore sizes, influenced by magnetic gradients, either at the pore or globally, this decay time indicates sand with low amounts of magnetic impurities or iron oxides at the grain surface.The first aquitard is imaged close to the expected depth as a conductive layer with less water content and decreased decay time, but the thickness appears too small.However, the error bar is so big that the layer can hardly be detected.The second aquifer exhibits increased resistivity compared to the aquitard but still lower resistivity compared to the first aquifer and a slightly lower water content.The latter might be explained by either a higher degree of compaction or a compensation effect due to a too large thickness.
Remarkably, the decay time of the second fine sand aquifer is far lower than for the first, even though the lithology is similar according to the borehole description.Even though the description of the latter states "weakly silty", an increased amount of small grains could decrease the medium pore size seen by the protons.This is partly supported by higher gamma ray intensities (see Fig. 2b), but decreased T * 2 due to increased amounts of fine-grained material cannot be proved without samples.Note that the T * 2 decay time is not only controlled by pore size but also by magnetic field inhomogeneities.Thus, faster decay may also be caused by magnetic gradients at pore scale due to magnetic impurities or iron oxides at the grain surface or in the fluid.Grunewald and Knight (2011) showed this influence as decreasing T * 2 with increasing magnetic susceptibility of the grains.Such an increase in magnetic susceptibility is partly supported by the brownish colour, but cannot be proved unless samples are measured.Nevertheless, since the genesis of this layer is different from the primary aquifer, internal pore scale gradients due to magnetic impurities only in the second aquifer are possible.The sensitivity of the MRS-measured T * 2 decay time to magnetic gradients is a principal disadvantage compared to laboratory T 2 or T 1 decay times.Measuring T 1 , the longitudinal relaxation, will avoid this ambiguity once appropriate measuring and inversion schemes are established.
The last layer is very conductive and shows again decay times typical for sand.It comprises both the clay-sand interlayers and the underlying saltwater, which can not be further distinguished due to the bad resolution below 50 m, the good conducting zone.
The above comparison of water content, decay time and resistivity with ground truth shows a main advantage of the combined application of MRS and VES.Both T * 2 decay times and resistivity interpreted independently are ambiguous.Fast decay may be due to smaller pores or magnetic gradients.Low resistivity may be saltwater as pore fluid or increased clay content.However, low resistivity with high decay times uniquely indicates an effect due to saltwater, whereas small decay times with high resistivity indicate magnetic gradients.
Before the obtained quantities can be further used, we are interested in their accuracy.As described above, we vary each model parameter individually until the model response exceeds the error model around the response from the best model.The range of valid models is expressed by the error bars (see Fig. 3), which are relatively large for water content.The variation of the decay time is far lower and increases with depth, nevertheless the layer parameters corresponding to the geological units are significantly different.The best parameter is resistivity, which is resolved in tight bounds except for the first aquitard.
Generally, the NMR parameters are poorly determined for the unsaturated zone due to the small signal.Also, the thin aquitard has a bad resolution and can not significantly be distinguished from the neighbouring layers.The thickness values have a quite small uncertainty because here all three parameters combine their resolution.Only the last layer's boundary is not well defined due to the high conductivity, as well as the NMR parameters in the conductor.In total, the uncertainty analysis shows that for this medium quality data, as in this case, a quantitative analysis needs to be done with caution.

Calibration sounding OD33
The next two soundings are located at the southern boundary of the dunes, where the thickness of the freshwater lens is still large.We used the identical experimental setup (cf.Table 2).Since the results are very similar, we show only OD33, a sounding that was made next to a well (P-OD33) where pumping tests have been carried out (Sulzbacher et al., 2012).See section 4.1 for the calibration of hydraulic conductivity with NMR parameters from this sounding.We decided to use a five-layer case for inversion.In contrast to CL2, the data quality was so excellent (cf.Table 2) that the intermediate layer between the aquifer was needed to reach a χ 2 of 1. Data fit and results are shown in Fig. 4.
Due to the small distance to CL2, the lithology is expected to be similar.Accordingly, the first aquifer has almost identical properties (≈ 31 % porosity, 200 ms decay time and about 80 m resistivity).Below, there is also a shallow, conductive layer with decreased decay time, which is therefore interpreted as silt.In contrast to CL2, the decay times of the second aquifer show again values of about 200 ms instead of 70 ms.Obviously, there must be a difference either in the amount of fine grains or magnetic ions in the aquifer, e.g.due to the pumping activities close to CL2.More probable are lateral variations of the deposition regime, which need to be further investigated.
Looking at the uncertainties the overall parameter resolution is significantly better than for CL2.Exceptions are the NMR parameters of the unsaturated zone and the silt layer.Also, the upper boundary of the silt is fairly uncertain, whereas its thickness is much better determined.While the models at CL2 underestimated the thickness of this first aquitard, the position and thickness at OD33 agree better with the borehole information.Due to the improved (by a factor of two) MRS data quality and the large contrast in the decay times, the silt layer can be nicely resolved and supports resistivity in the joint inversion.

Prediction sounding SKD
The last sounding is located at the easternmost boundary of the freshwater lens close to the North Sea.In this flat area storm surges, as known from the storm Kyrill in 2007, are regularly adding saltwater from the top.In this vulnerable zone between deep saltwater intrusion and surface salinity, the dynamics depends highly on the distribution of the hydraulic conductivity.We used a smaller loop (25 × 25 m), shorter pulses and a higher number of pulse moments, resulting in a higher resolution for the shallow depths.After processing, a medium stacking error (Müller-Petke et al., 2011) of about 9 nV was determined.Corresponding to the maximum amplitudes of the data cube, this is a signal-to-noise of 32 (Table 2).Noise-cancellation did not improve S/N (signal-to-noise ratio) and was therefore not used.
Resistivity information could not be derived from a nearby VES site (cf.Fig. 1).However, by comparing two soundings at the same distance to the dunes we found only little differences pointing to 2-D conditions.Therefore, we used the closest measurement for the joint inversion, about 150 m south of the MRS coil.From neighbouring boreholes a shallow and thick silt layer was presumed, but we had no reliable idea about the subsurface layering.Therefore, we first conducted independent inversions of both data sets and increased the number of layers subsequently.MRS data could already be fitted to a χ 2 level of 1.2 using a three-layer case (see Fig. 5a,b).
Whereas water content is almost constant at about 30 %, the decay times showed a significant decrease from 200 to 75 ms in the second layer, which can be interpreted as the silt layer.On the other hand, four layers (Fig. 5c) were needed to fit VES data, which is obvious from the apparent resistivity curve (Fig. 5f).The actual values are very low, hinting to high salt concentrations in the fluid.Interpretation of the resistivity model alone, however, is hard since silty or clayey layers have similar resistivity ranges as salt or brackish water in sand.
The position of the second layer boundary is amazingly equal, which gives hope for a successful joint inversion.Again, we started with a small number of layers, but a significant decrease of the fit close to the error level was only possible with 5 layers (χ 2 = 1.07), as already indicated by the independent models.Figure 6 shows the joint inversion result.The model responses can hardly be distinguished from the single inversions.
A sandy layer (200 ms decay time) is again on top with already low resistivity originating from the flooding.The good conductor below is, as in Fig. 5, split into two different decay times, most probably a sandy layer with 200 ms and a silt layer with 75 ms.The silt obtains slightly higher water content, which is very plausible, although the uncertainties bigger than the overall variations.The uncertainties of the second and third layers signal that they cannot be distinguished with resistivity methods, but likely can be by the T * 2 time.
Below this silt layer there is fine sand again with large decay times, but a resistivity jump from about 20 down to 2 m.Whereas the lower represents saltwater, we interpret the upper as laterally influenced by the freshwater lens, which is much shallower than in the dunes and more conductive due to lateral diffusion of saltwater.Looking at the uncertainties we see that the best determined parameter is decay time, followed by resistivity.Water content, although nicely detected by MRS, cannot be used for discrimination.The uncertainties of the layer thicknesses increase with depth but are very small due to the fact that all three parameters contribute to them.In total the model can be considered very reliable even if no ground truth can be used.Several facts are responsible for this: the very low absolute noise level and the large amount of stacks, the number of pulse moments, and the short dead time that makes a clear detection of fastdecaying water possible.

Derivation of hydraulic properties
The Borkum groundwater model by Sulzbacher et al. (2012) requires the distribution of the following input parameters: porosity, current salt concentration and hydraulic conductivity.All three obtained parameters can be derived from the quantities obtained in the joint inversion.
Porosity correlates with the NMR water content in case of full saturation.However, the detected water must be interpreted as mobile water, since adhesive water in very small pores leads to an underestimation of the total porosity if clayey material is present (Lubczynski and Roy, 2005).Vouillamoz et al. (2007) introduced an additional calibration factor for MRS water content to total porosity.This calibration is derived from pumping tests via specific yield in the case of an unconfined sandy aquifer.We do not expect large amounts of clayey material, and thus we avoid another calibration factor.
Salt concentration is deduced from a modified Archie equation (Waxman and Smits, 1968), relating fluid conductivity σ f and bulk conductivity σ b : where σ s is a surface conductivity originating from finegrained material and F is the formation factor, which is related to porosity by the cementation factor m.Even though there are widely used literature values, σ s and m can be obtained on-site.
The data of about 15 direct push soundings in the survey area were used.These soundings measured both bulk resistivity by a four-point measurement and took fluid samples for retrieving fluid conductivity.We used 25 data pairs from the first aquifer showing a wide range of salt concentrations, and curve-fitting resulted in σ s = 3.66 mS m −1 and F = 3.83.From the obtained porosity of the first aquifer, the latter is equivalent to a cementation exponent of m=1.26, very close to literature values of 1.3 for loose sand.These values are subsequently used to derive the fluid conductivity from porosity φ and resistivity ρ and further deduce a NaCl concentration.This procedure is in contrast to Vouillamoz et al. (2007), who used a rough linear equation with a single calibration factor to estimate the conductivity of the pore fluid from bulk conductivity due to expected clay contents that prohibit the use of Archie's law.
Hydraulic conductivity K is obtained from porosity φ and decay time T * 2 using a semi-empiric equation from Kenyon (1997):  tortuosity but also NMR rock parameters such as surface relaxivity.We use an early proposed set of factors from Seevers (1966), validated on quartz powder and sandstones (small porosities), both measured in laboratory conditions (T 2 instead of T * 2 ).However, C S is still a calibration factor but different from Eq. ( 6).Several authors have used this equation and observed appropriate values for C S between 30 × 10 −4 and 326 × 10 −4 m s −3 (Mohnke and Yaramanci, 2008).
Aside from grain size analyses or flow experiments, a standard method at the aquifer scale is to carry out pumping tests.Several of these were conducted in the east of Borkum in water test wells.Additionally, fluid conductivities are wellknown.We chose the well P-OD33 for calibration since it represents a simple situation with a typical fine sand aquifer that is far away from pumping wells.
For the upper aquifer at P-OD33, a transmissivity of 9.86 m 2 s −1 was determined from a pumping test (Sulzbacher et al., 2012).Taking the determined thickness of 14 m into account, this corresponds to a hydraulic conductivity of K = 7.04 × 10 −5 m s −1 .By inserting the porosity φ = 32.3 % and decay time T 2 = 215 ms into Eq.( 7), we obtain a calibration factor of C S = 47×10 −4 m s −3 , which is well within the literature range (Mohnke and Yaramanci, 2008).
Since OD33 was already used for calibration of hydraulic conductivity and CL2 exhibits a poorer data quality, we show the derivation of the hydraulic parameters for the sounding SKD.As described, each two of the three primary parameters are combined to obtain hydraulic conductivity and fluid conductivity.The latter is expressed as total dissolved solids (TDS) using a Chloride (Cl − ) concentration conversion factor derived by Sulzbacher et al. (2012) from water sample analyses.We did not apply the Archie equation for the silt layer, since we do not know its surface conductivity.The results are summarized in Table 3.
The hydraulic conductivities of the fine sand layers are in a very plausible range between 4 × 10 −5 and 7 × 10 −5 m s −1 .
Only the last layer exhibits a high value due to the unusually high T * 2 .The TDS concentrations of the lowest layer are close to that of seawater.The same holds for the second layer, where saltwater inserted by flooding events is obviously assembling on top of the silt which is acting as a semi-permeable barrier.Both the top layer and the layer beneath the silt show brackish conditions.

Uncertainty of the derived parameters
As the primary results of any physical experiment require a measure of reliability or uncertainty, the final outcome of our survey requires uncertainties of the determined hydraulic parameters φ, K and TDS.We apply the maximum error propagation, i.e. the uncertainties caused by all individual values add up to the error of the target value.
The retrieved TDS concentrations obtain the same relative error as the fluid conductivity or resistivity.Furthermore, we assume the concentration exponent and surface conductivity are constant.The first should not vary that much and the latter plays only a minor role in sandy sediments.From Eq. ( 5) we can directly derive by differentiation As a first order approximation, the relative error is slightly higher than the sum of the relative errors of both primary parameters.Similarly, we derive from Eq. ( 7) Of course the parameters are not completely independent, e.g. an overestimation of water content is often accompanied by an underestimation of decay time.However, we can use the result as a conservative guess.The first factor can only be retrieved from the calibration itself and is of the same order of the rest (relative error in porosity and twice the relative error in decay time).If several calibration wells are available, the uncertainty in C S can be drastically decreased and, additionally, be estimated from the variation of the retrieved values.
We apply Eqs. ( 8) and ( 9) to the results of SKD and assume the calibration factor known in order to show what we can achieve in case of a good sounding.The relative deviations of all parameters are summarized in Table 4.
The relative primary parameter variations are all between half an order of magnitude, most of them at about 10-20 %.The worst values are obtained for the silt and the lowermost layer, where the resolution is low.Consequently, for these also the secondary parameters are not well resolved.However, the three sandy layers show very small ranges due to the good data quality.Vouillamoz et al. (2007) presented slightly but generally higher uncertainties.There are some reasons for this.(i) As  2007) calculated the uncertainty from statistics of many field MRS and VES measurements in contrast to the uncertainty analysis of a single sounding.

Conclusions and outlook
We have presented a new methodology for investigating layered groundwater systems.It combines the analysis of a vertical electric sounding and a magnetic resonance sounding by using a block model.The retrieved parameters resistivity, water content, and decay time allow for determining much less ambiguous lithological models compared to resistivity soundings alone, since the decay time is sensitive to pore size and not to fluid conductivity.The technique is applied to measurements on the North Sea island of Borkum, a freshwater/saltwater system possibly threatened by ongoing climate change, and yields a clear lithological description that was not available from noninvasive methods so far due to the ambiguity of their results.The geometric coupling of the two independent methods leads to reduced uncertainty of all parameters and layer thicknesses as compared to single inversions.
The methodology can easily be modified using any other resistivity sounding method, such as time domain or frequency domain electromagnetics, or any combination of them.We expect further improvement once the longitudinal relaxation time T 1 can be reliably measured and reconstructed since T 1 is not influenced by magnetic gradients in the aquifer.Enhancing the method to 2-D or 3-D, either by 2-D/3-D surveys or lateral/spatial constrained 1D data, would be necessary to describe whole catchments.Alternatively, available airborne EM data and geostatistical analyses could be used for generating 3-D models.
From the obtained parameters we are able to derive the parameters porosity, hydraulic conductivity, and salt concentration, which are needed for a density-driven modelling of the hydraulic system.Key issues are the used petrophysical relations and their coefficients.The conversion from bulk conductivity to fluid conductivity is based on the analysis of direct-push measurements.The hydraulic conductivity is obtained by a semi-empiric relation after calibrating with a pumping test.Although calibration is still necessary and bears additional uncertainty, the multi-parameter process is superior to deriving TDS concentrations from electromagnetic or DC resistivity measurements alone.
If the methodology is further developed and carried out more routinely, it will be able to improve the hydraulic models both by a more accurate geometrical description and more reliable estimates of hydraulic conductivity and current salt concentrations.As a result it will significantly enhance the forecast of climate changes on groundwater systems and thus contribute to sustainable water management.

Fig. 1 .
Fig. 1.Measuring area in the east of Borkum and locations of MRS measurement (red) and reference (black) coils, VES midpoints (blue tri-stars), the CLIWAT II borehole (magenta) and HEM soundings (red dots) in the east dunes of Borkum.

Fig. 4 .
Fig. 4. Joint inversion results (subplots as in Fig. 3) of the MRS and VES data for the location OD33.

Fig. 5 .Fig. 6 .
Fig.5.Results of independent inversions (subplots as in Fig.3) of the MRS and VES data for the location SKD using 3 and 4 layers, respectively.

Table 1 .
Lower and upper parameter bounds for the individual parameters.

Table 3 .
Retrieved hydraulic parameters for sounding SKD.

Table 4 .
Vouillamoz et al. (2007)all primary and secondary parameters from sounding SKD.Vouillamoz et al. (2007), including a fixed geometry, here derived from MRS, decreases uncertainty.Our joint inversion of MRS and VES combines the best of both methods and improves resolution, and therefore decreases uncertainties.(ii)QTinversion as shown by Mueller-Petke and Yaramanci (2010) is a more general and accurate approach.The block-QT inversion reduces the number of unknowns and therefore the accuracy of its estimation.(iii)Noise cancellation techniques (Müller-Petke and Yaramanci, 2011) improved the data quality of MRS significantly.The data presented have a very good signal-to-noise ratio.(iv)Vouillamozet al. (