A further development of a special inversion method for the magnetic resonance sounding method is presented. It allows noninvasive parameter estimation for water retention models by parameterizing the capillary fringe if the water table position is known and included in the inversion process.


The magnetic resonance sounding (MRS) method is usually applied for delineation and characterization of aquifer system stratification. Its unique property, distinct from other hydrogeophysical methods, is the direct sensitivity to water content in the subsurface. The inversion of MRS data yields the subsurface water content distribution without need of a petrophysical model. Recent developments in instrumentation, i.e., decreased instrumental dead times and advanced noise cancellation strategies, enable the use of this method for investigating the vadose zone. A possible way to interpret MRS measurements with focus on water retention (WR) parameters is an inversion approach that directly provides WR parameters by modeling the capillary fringe (CF inversion). We have developed this kind of inversion further to account for different WR models and present a sensitivity study based on both synthetic and real field data. To assess the general applicability of the CF inversion, we analyzed the resolution properties for different measurement layouts and the parameter uncertainties for different realistic scenarios. Under moderate noise conditions and if the water table position is known, all WR parameters except the residual water content can be reliably estimated. The relative accuracy of the estimated pore distribution index estimation is better for larger CF. Small measurement loops of 5-m diameter achieve the best resolution for shallow investigation depths of <10 m.

Besides laboratory- or point-scale methods, there are few nondestructive geophysical methods sensitive to moisture distribution in the vadose zone such as ground-penetrating radar (GPR) in different setups (Huisman et al., 2003). For example, GPR reflections from the water table are sensible to water retention parameters (Saintenoy and Hopmans, 2011). Those parameters can also be estimated in a coupled hydrogeophysical inversion using time-lapse GPR data (e.g., Scholer et al., 2013). Another method that can provide water retention parameters within a time-lapse framework is electrical resistivity tomography (e.g., Mboh et al., 2012).

The surface nuclear magnetic resonance (NMR) method is commonly used for groundwater investigation and aquifer characterization. Its unique and beneficial property, distinct from most other hydrogeophysical methods, is the direct sensitivity to the subsurface water content (Legchenko et al., 2004; Yaramanci and Müller-Petke, 2009; Hertrich et al., 2009). Petrophysical models for estimating the water content from geophysical parameters such as electrical resistivity or dielectric permittivity exhibit additional error sources and lead to higher uncertainty of the final results.

The one-dimensional application of surface NMR is often called magnetic resonance sounding (MRS). In addition to water content, the method can provide estimates of hydraulic conductivity if adequate calibration data from pumping tests (Legchenko et al., 2004) or core sampling is available (Mohnke and Yaramanci, 2008). First attempts to investigate the unsaturated zone were made about 10 yr ago (Lubczynski and Roy, 2003; Roy and Lubczynski; 2005). They showed that pore water above the water table in an unconfined aquifer can be measured with MRS and can be interpreted qualitatively. However, the MRS instruments of the first generation were not optimized for measurements in the shallow subsurface, i.e., at a depth of <10 m. Moreover, the applicability of MRS is often limited by a low signal/noise ratio due to its very small signal strengths of <100 nV for measurements in partially saturated material. Therefore, the potential of the method for investigating and characterizing the unsaturated zone could not yet be exploited.

A new generation of devices is now available that allows improved noise reduction techniques on the one hand (Walsh, 2008; Dlugosch et al., 2011; Müller-Petke and Costabel, 2013) and an optimization of the measurement layout with focus on the shallow subsurface (Costabel and Yaramanci, 2011a; Walsh et al., 2011) on the other hand. The technical development of the MRS equipment to particularly improve its applicability for investigating the unsaturated zone is an ongoing process (Walsh et al., 2014). Costabel and Yaramanci (2011a) suggested an approach for observing the water content distribution in the capillary fringe (CF) of unconfined aquifers. This approach is based on the assumption that the CF is in hydrostatic equilibrium and can be approximated using water retention (WR) functions. They used the WR model of Brooks and Corey (1964) and showed that the MRS data inversion can be modified such that WR parameters are directly provided. In this first study, however, only one WR model was investigated and only one special loop configuration was considered. In this study, we developed this approach further with an improved inversion scheme that also allows the use of the models of van Genuchten (1980) and Kosugi (1996). We investigated and assessed the applicability of this approach with a detailed error analysis based on simulated and real field measurements. We made use of different measurement layouts and took into account different water retention models.

Basics of Magnetic Resonance Sounding

The MRS method is based on the NMR phenomenon, i.e., the nuclear magnetic excitation and observation of hydrogen nuclei 1H in water molecules (e.g., Levitt, 2002; Legchenko and Valla, 2002). A 1H nucleus has a magnetic dipole moment caused by its spin angular momentum, which is aligned with the local predominant static magnetic field B0. In the case of MRS, this is the Earth’s magnetic field. For conducting MRS, a circular or square cable loop is placed at the surface with a diameter or side length from a few meters for shallow investigations up to 150 m for deeper aquifers. During the transmission mode, the loop is energized with an electromagnetic pulse, i.e., for a short period of time (5–40 ms) an alternating current is produced. The product of the current strength and the pulse length is referred to as the pulse moment, q. The pulse causes an electromagnetic (EM) field that propagates into the subsurface. The frequency of the energizing EM field is chosen such that it equals the precessing frequency of the nuclear spins, i.e., the Larmor frequency (Levitt, 2002). Consequently, the spins are forced to tilt away from their equilibrium positions. After the pulse is switched off, the spins relax back to their state of equilibrium, which induces a measurable voltage in the surface loop alternating with the Larmor frequency and decaying at a certain relaxation rate. The voltage is recorded after a device-specific dead time, which is required to switch the device from the transmission mode to the receiving mode. The recording can be made using either the transmission loop or a second receiving loop. The envelope of the recorded signal is, in general, a multiexponentially decaying signal: 
where t is the time after the pulse, T2* is the relaxation time, and I is the intensity function corresponding to T2*. The initial amplitude 
of the recorded signal is proportional to the number of excited 1H nuclei, i.e., to the water content θMRS in the sensitivity range of the MRS measurement: 

The subscript MRS indicates that we consider only the water content measurable with MRS. Generally, the NMR relaxation time increases with pore size (Coates et al., 1999) and saturation degree (Ioannidis et al., 2006). Consequently, NMR signals from water in small pores and residual water in drained large pores with T2* times shorter than the instrumental dead time are undetectable (Boucher et al., 2011; Knight et al., 2012). This means that capillary and clay-bound water is normally not measured with MRS. Actually, the effective dead time is even longer than the instrumental dead time because the relaxation process already starts during pulse application (Walbrecker et al., 2009). When determining θMRS from MRS measurements, the effective dead time is accounted for by taking the middle of the pulse as t = 0 (Walbrecker et al., 2009). After this correction, E0 represents the initial amplitude, and commonly, the corresponding θMRS is considered to be the effective (free) water content (Legchenko et al., 2004). In practice, however, it is expedient to verify this assumption for the investigated material at the site of interest.

The NMR devices for application in boreholes or in the laboratory are designed such that they can observe the relaxation process with high accuracy either parallel (longitudinal or T1 relaxation) or perpendicular to B0 (transverse or T2 relaxation). The corresponding relaxation times are called T1 and T2, respectively. Both parameters are strongly correlated to the pore size of porous media (Coates et al., 1999). Characterizing the multiexponential character of T1 or T2 signals from saturated porous media allows the reconstruction of pore size distributions for sandstones and limestones (e.g., Kenyon, 1997) and water retention curves for sandy soils (Costabel and Yaramanci, 2013). Also, the T2* relaxation measured by MRS depends on the pore size (Mohnke and Yaramanci, 2008). Unlike the T2 relaxation, however, it is additionally affected by inhomogeneities in B0 (Grunewald and Knight, 2011). Therefore, and because of the lower data quality, T2* measurements will probably not allow the estimation of pore size distributions in many cases.

When formulating MRS signals in a multiexponential form (Eq. [1]), we account for the fact that one single MRS measurement covers a certain depth range that can, in principle, consist of more than one relaxation regime, e.g., different aquifers with different mean pore sizes. The resulting signal is then a superposition of the signals from all layers in the covered depth range. In other words, in addition to θMRS, T2* is a function of depth z as well (Müller-Petke and Yaramanci, 2010): 

The function κ(q,z) is the sensitivity function of the MRS, i.e., it describes the observed depth ranges of the MRS measurements. In practice, κ(q,z) is calculated by modeling the EM field of the loop, for which the subsurface distribution of the electric resistivity must be known (Legchenko and Valla, 2002; Valla and Legchenko, 2002). The depth of investigation for a single MRS measurement is controlled by q. An increasing q shifts the observation depth to deeper regions. The entire MRS experiment is an ensemble of single measurements with increasing q. During the measurements, the increase in q is usually realized by increasing the transmitting current while the pulse length remains constant. Finally, the interpretation of an MRS data ensemble is done by solving Eq. [4] for θMRS(z) and T2*(z), i.e., the distributions of water content and relaxation time with depth (Müller-Petke and Yaramanci, 2010). Due to its ill-posed nature, the inversion must be regularized. This can be done by introducing smoothness constraints (Müller-Petke and Yaramanci, 2010) or a fixed number of layers for a so-called block inversion (Günther and Müller-Petke, 2012).

When the investigation is focused on the distribution of θMRS(z) only, Eq. [4] can be simplified to 

In such cases, the inversion problem is solved in two steps (Legchenko and Shushakov, 1998). First, the single measurement curves are fitted exponentially to get the initial amplitudes E0. Second, Eq. [5] is solved by inversion to identify θMRS(z).

Magnetic Resonance Sounding for Vadose Zone Investigations

A Field Example

Figure 1a shows the data from an MRS example conducted in April 2013 at the test field Barnewitz/Nauen west of Berlin, Germany. Beneath the topsoil, a homogenous sand layer down to a depth of approximately 7 m is present (Costabel and Yaramanci, 2011a). This layer represents an unconfined aquifer. The water table ztable = 1.2 m was measured in a borehole at a distance of 20 m from the position of the MRS loop. The MRS measurement layout was focused on the vadose zone: a small circular measurement loop with a diameter of 6 m was used (transmitter loop with two turns and receiver loop with 24 turns). In addition to the measurement loop, two additional loops were placed at a distance of 30 m to measure the x and y components of the surrounding electromagnetic noise. In the post-processing, the data from these loops were used to reduce the noise level from the MRS measurements by harmonic noise cancellation (Müller-Petke and Costabel, 2013). The MRS data set consists of 48 pulse moments ranging from 0.006 to 0.21 A s. The instrumental dead time was 2.5 ms and the pulse length 9.5 ms yielding an effective dead time of 7.25 ms.

Three different inversion results from this data set, after inverting Eq. [4] with smoothness constraints, are given in Fig. 1b and 1c, showing the distributions of θMRS(z) and T2*(z), respectively. Each inversion was performed with an individual degree of smoothness for the resulting model, which is controlled by the so-called regularization parameter λ. A larger λ enforces higher smoothness in the model at the expense of a worse data fit (Aster et al., 2005). To decide whether the data fit is acceptable or not, we used the χ2 value (error-weighted misfit, after Günther and Müller-Petke, 2012) that should be close to 1 for correct error estimates. For the inversion, the choice of λ must be made with care; choosing its value too small often causes artifacts in the resulting model. This effect is demonstrated by the green curve in Fig. 1b. In the corresponding model, the topsoil appears to have a θMRS > 0.3, which supports nearly saturated conditions. However, this is very unlikely for a sandy soil. After performing the MRS measurements, we took soil samples down to a depth of 50 cm covering the topsoil layer to verify the MRS inversion result. We found that the actual θ of the topsoil was in the range of about 0.20 to 0.22 (symbols in Fig. 1b). The black curve (λ = 50,000) in Fig. 1b seems to be a reliable model, although the fit to the data is slightly worse than for the green curve (λ = 10,000). When λ was further increased, the resulting model (blue curve in Fig. 1b) appeared to be overregularized, i.e., the smoothness was too strong and the data fit worse. The described principle is also true when considering the resulting T2*(z) distributions depending on λ in Fig. 1c. For our field example, however, we did not have reference values for T2*, so we could not decide which model in Fig. 1c is the most reliable one. The green curve seems to be underregularized (like the green curve for θ in Fig. 1b) because changes in T2* beneath ztable in the homogenous sand layer are very unlikely. Apart from that, T2* generally increased with θMRS, which is in line with expectations.

For the approach investigated in this study, we neglected the NMR relaxation and used only the initial amplitudes of the MRS measurements for further analysis. In other words, we concentrated on the interpretation of the water content distribution only. Costabel and Yaramanci (2011a) suggested constraining the MRS inversion by describing the θMRS increase in the capillary fringe with WR models. We refer to this approach as the capillary fringe (CF) inversion. For its application, only the transition between the unsaturated and the saturated zone is considered, i.e., the topsoil is neglected. To apply the CF inversion to the MRS data set presented in Fig. 1a, the NMR signals corresponding to the topsoil, i.e., the smallest q must be excluded. Figure 1d shows the same data as in Fig. 1a, but now without the signals for q < 0.9A ms. The inversion runs were repeated with the same λ values as described above and the results are depicted in Fig. 1e and 1f. We note again the varying smoothness of the resulting models, but the increased water content in the topsoil does not appear, even for the lowest λ. All models show in principle the same trend: the θMRS(z) distribution starts with the smallest values in the vicinity of the subsurface and increases with depth until saturation, as expected for the CF. However, none of these models can reliably reconstruct the CF because all of them overestimate the depth of the saturated zone. We conclude that the usual smooth inversion scheme is not suitable for characterizing the CF properly.

Reformulation of Forward Model and Capillary Fringe Inversion

For a plausible physical description of the CF, we used common WR models. Water retention models describe the volumetric water content θ as a function of the capillary pressure head h. The pore space of soils is assumed to behave like a bundle of capillaries and thus WR functions represent cumulative capillary radius distributions. With decreasing radius, the absolute value of h increases and vice versa. In contrast to the hydrostatic pressure head below ztable, the capillary pressure head above ztable is defined with a negative sign to account for the direction of the capillary force against gravity. In equilibrium, the absolute value of h equals the height above ztable, which is characterized by h = 0. For the CF inversion, we consider the CF in equilibrium. Consequently, θ inside the CF is determined by the chosen WR model and decreases accordingly with increasing distance |h| to ztable. For the reformulation of the MRS forward problem, we modified the general formulation in Eq. [5] by substituting θMRS(z), with θ described by WR functions. In doing so, we accepted that the modified MRS forward problem will describe solely the transition between the unsaturated and saturated zones. In practice, one must therefore carefully determine the MRS measurement conditions to assure that the whole CF is covered adequately by the measurements (Costabel and Yaramanci, 2011a). On the other hand, the assumption of a homogenous, unconfined aquifer must be plausible.

Regarding the choice of WR models for our approach, we decided to reduce the degree of ambiguity for the parameter estimation as much as possible, i.e., we chose models with as few parameters as possible. In addition to the model of Brooks and Corey (1964), hereafter referred to as the BC model, we took the models of van Genuchten (1980) and Kosugi (1996) into account, from now on abbreviated with VG and KO, respectively. Each of these models consists of four parameters (Table 1): θS and θR (both dimensionless) are the volumetric water content at saturation and the residual water content, respectively; the other two parameters describe the transition between θS and θR, and their impact is similar in each model. The parameter h0 shifts the transition in relation to ztable (h = 0), the other one defines the slope of the transition. Regarding its physical meaning in each model, we refer to the latter as the pore distribution index (PDI) for simplification. However, the exact meaning of h0 and the PDI is slightly different in the individual model. For the BC model, h0 represents the air-entry point, which corresponds to the height of the saturated zone, i.e., the saturated region above ztable. The PDI λ appears to be a simple exponent in the BC model. For the VG model, h0 and n represent purely empirical parameters. The KO model describes a lognormal distribution function, with h0 representing the mean value and σ the standard deviation.

The WR models are defined by distinguishing between two cases, h ≥ 0 and h < 0. In Table 1, this distinction is expressed in relation to ztable, where h = 0, as it is meaningful for the MRS interpretation. It is obvious that ztable must, in principle, be considered as an additional free parameter to be fitted during the inversion. However, ztable and h0 cannot be estimated simultaneously because MRS is sensitive to θ but not to water tension. Thus, it can only provide the depth at which the saturated zone begins (e.g., at ztable–h0 for the BC model). As MRS provides no indication of the position h = 0, any corresponding attempt to fit both ztable and h0 at the same time would lead to a high degree of ambiguity and would therefore destabilize the numerical optimization process. Therefore, one of the two quantities must be given a priori and must be fixed during the inversion.

The inversion was realized by the GIMLi library of Günther and Rücker (2009) and was based on a Levenberg–Marquardt optimization algorithm. The WR parameters were kept within reliable boundaries using a cotangens model parameter transformation (Dlugosch et al., 2013). Boundary and starting values for all inversion runs are given in Table 2. We set θR = 0 for all inversion runs in this study because θMRS is commonly considered to be the free water content, i.e., the water content reduced by θR. Nevertheless, we formulated the forward models in Table 1 considering a non-zero θR to define a general basis for the CF inversion. As explained above, the assumption θMRS = θ − θR cannot be considered to be the general case and should be evaluated at the site of investigation. Moreover, future development of MRS technology may allow measurements with further shortened dead time and may enable the estimation of θR as well.

Study with Synthetic Data

Simulation of Realistic Measurements

To investigate and assess the potential of MRS for estimating the WR model parameters, we simulated realistic MRS measurements numerically. Two subsurface water content models were considered, one representing a sandy soil (Model 1) and the other a loamy sand (Model 2). For both models, θR and θS were set to 0 and 0.35, respectively. The increasing θ(z) of the capillary fringe was determined by VG parameters according to the database of Hammecker et al. (2004): h0 was 20 cm for both models, n = 2.3 for Model 1, and n = 1.7 for Model 2. Figure 2a shows the corresponding θ(z) distributions for the case of ztable = 3 m. In Fig. 2b, the κ(q,z) functions (i.e., sensitivity functions, see Eq. [4]) of the simulated MRS are shown. These functions were calculated for a loop with a 10-m diameter and with 40 pulse moments ranging from q = 0.01 to 1 A s. With the numerical simulations, we wanted to test and compare the performance of the CF inversion with regard to different situations. The study included different noise scenarios and measurement layouts (varying number of pulse moments nq and loop sizes), as well as variations of ztable. For all simulations in this study, we chose the q configuration such that the main part of the sensitivity function corresponding to the highest q was inside the saturated zone, i.e., beneath ztable, whereas the smaller pulse moments covered the unsaturated zone with minimum q = 0.01 A s. Figure 2c shows the simulated E0(q) according to Eq. [5] for the θ(z) models in Fig. 2a, contaminated with Gaussian noise of 1-nV standard deviation. It should be noted that this noise level refers to the uncertainty in E0 in real data sets, i.e., the error in E0 after the exponential fitting of the signals. It is not identical to the actual mean noise amplitude on the recorded time series.

As described above, the CF inversion cannot reliably be applied when approximating both ztable and h0 simultaneously. At least one of the two must be pre-estimated and held fixed during the inversion. For the CF inversion result shown in Fig. 3, h0 was given a priori. In this figure, we show the benefit of the CF inversion compared with a simple block inversion that considered the saturated and the unsaturated zones as simple layers with a sharp interface. Figure 3a shows the simulated data of Model 1 as in Fig. 2c, together with the approximations. The CF inversion was based on the VG model. For both inversions, the root means square error (RMS) was similar to the noise level of 1 nV, which means that both lead to equivalent results within the data error. In Fig. 3b, the inversion results from the block and CF inversions are shown together with the original model. Both models estimated θS correctly, but obviously the ztable estimation of the block inversion was underestimated because the transition between the saturated and unsaturated zones remained unconsidered. The CF inversion reconstructed the model accurately, and the ztable estimate is correct. Figure 3c and 3d show the objective functions considering the parameters ztable and θS for both inversion processes. For both, the minimum occurs without additional local minima, and the size of the minimum is similar, indicating a stable numerical process and a similar scattering of the parameter estimates. However, the minimum for the block inversion is shifted towards smaller ztable. Consequently, a systematic underestimation must occur. When using the CF inversion, this underestimation can be avoided. In addition, WR parameters can be estimated, which is demonstrated and assessed below.

Accuracy of Water Retention Parameter Estimates

From now on we assume that ztable is given as in the real case studies below. If this is the case, the CF inversion estimates the WR parameters, i.e., θS, h0, and the PDI. In principle, also θR could be estimated. However, throughout this study we manually set this parameter to zero and neglected its influence. So far, all of our real measurements have indicated that the residual water content cannot be measured with MRS or is systematically underestimated (Costabel and Yaramanci, 2011a, 2011b). Figure 4 shows a simulated MRS (Fig. 4a) based on Model 1 and the corresponding CF inversion results (Fig. 4b) using all three considered WR models from Table 1. Figure 5 shows the same for Model 2. For these simulations, the same measurement parameters as for Fig. 3 were used, i.e., a 10-m loop, 40 pulse moments in the range from 0.01 to 1 A s, and 1-nV noise. Above, Model 1 and 2 were defined using the VG parameterization. We also used the BC and the KO models. Thus, we redefined the parameterization of Model 1 and 2 with regard to BC and KO. In doing so, the BC and KO parameters were determined manually such that the resulting WR curve approximated the original VG model best. For BC, the parameters were λ = 1.4 and h0 = 20 cm (Model 1) and λ = 0.7 and h0 = 20 cm (Model 2). Regarding the KO model, we used σ = 1.0 and h0 = 30 cm (Model 1) and σ = 1.7 and h0 = 50 cm (Model 2).

Regarding the CF inversion for all three models, ztable was held fixed and the parameters θS, PDI, and h0 were estimated. As is obvious from Fig. 4 and 5, the inversions using the different WR models led to similar θMRS(z) distributions. To assess the accuracy of the parameter estimation, we analyzed their model covariance matrices. These are depicted in Fig. 4c to 4e (and Fig. 5c–5e for Model 2) after normalization to diagonal elements showing values between −1 (negative correlation) and 1 (positive correlation). We note that θS is weakly correlated with the other two parameters for each model, whereas the relationship between the PDI and h0 is characterized by a strong positive correlation for the VG and BC models and a strong negative correlation in the case of the KO model. This means that the θS values will be estimated with higher accuracy regardless of the estimation of the other two parameters, whereas an error in the PDI estimate will automatically lead to an error in the h0 estimate and vice versa. In addition to the correlation matrices, the square roots of the diagonal elements of the corresponding covariance matrix are plotted in Fig. 4c to 4e (and Fig. 5c–5e for Model 2). These values represent rough estimates of the uncertainties for the parameter estimates (Aster et al., 2005). They might be underestimated because the cross-dependencies between the parameters are not accounted for. An important objective of the following numerical study was the assessment of whether the diagonal elements are reliable uncertainty estimates. When comparing the results of the diagonal elements of Models 1 and 2, we note that the PDI for Model 2, i.e., the loamy sand model with a smoother retention curve, was estimated with higher accuracy for all models, whereas the uncertainties of the other parameters were similar. Note that the relative uncertainty of σ (i.e., the PDI of KO) was higher for Model 2 than Model 1. Because we did not observe significant differences in the performance of the three WR models, the following detailed analysis was done using the VG model only.

We next tested the influence of the noise level and the number of pulse moments nq on the CF inversion result using the same setting as above. The noise level was varied from 0.1 to 4 nV (x axes in Fig. 6 and 7). Three different numbers of pulse moments nq were used: 20, 40, and 80 (top, middle, and bottom of Fig. 6 and 7). For each combination of noise level and nq, we repeated the simulation 100 times and present in Fig. 6 and 7 the parameter estimates for Models 1 and 2, respectively. In doing so, we chose a visualization using the well-known box plot: each subplot in Fig. 6 and 7 shows the median (black line), the quartile (dark gray area), and the range where outliers occurred (light gray area) as functions of the noise level on the x axis. Additionally, the square root of the predicted variances, σpred, is plotted with dashed lines. We observed a general increase in uncertainties with increasing noise level and decreasing nq. For both models, θS was the parameter with the best performance. The scattering of the θS estimates ranged from less than ±0.01 for a low noise level up to approximately ±0.05 for the highest noise level. For both models, σpred is an effective and reliable measure of the real uncertainties. For the parameters h0 and PDI, we note that the scattering is asymmetric, i.e., outliers tend to overestimations. The asymmetric character of the distribution increases with the noise level. For Model 1, the h0 parameter increasingly tends to systematic underestimation with increasing noise level. In contrast, for Model 2, a systematic underestimation cannot be observed. Instead, the quartile for noise levels >0.5 nV is three times larger than for Model 1, and also the range of the occurring outliers is increased. This means that h0 was generally estimated with higher uncertainty for larger capillary fringes. A similar behavior was reported, e.g., by Schaap et al. (2001): with decreasing PDI n (corresponding to a larger CF in this study), they found an increasing uncertainty of parameter α (= reciprocal of h0 in the VG definition of this study). Regarding the uncertainty of the PDI, the opposite is the case, i.e., with increasing PDI its relative accuracy increases. For Model 1, the n estimates are systematically overestimated at noise levels >1 nV (for nq = 20 and 40) and >2 nV (for nq = 80). In contrast, for Model 2 the n estimates are reliable even for the higher noise levels, which is indicated by a constant median and a small increase in the quartile with increasing noise. However, the risk of outliers overestimating the n parameter increases with increasing noise level for both models, whereas underestimations of n are unlikely. We note that σpred represents the overestimations of n and h0 rather than the underestimations.

For the next test, we varied the loop size in three steps (5-, 10-, 20-m diameter; top, middle, and bottom rows in Fig. 8 and 9) and ztable (range 1–10 m, x axes in Fig. 8 and 9) and analyzed the corresponding influence of these parameters on the accuracy of the CF inversion. We expected different MRS loop sizes to have individual resolution properties at different depth ranges. For the simulations, the noise level and nq = 40 were fixed. Under realistic measurement conditions, receiver loops of different sizes will exhibit different absolute noise levels at the same location. The larger the loop area, the higher is the magnitude of the induced noise. To compare the simulated inversion results from different loop sizes adequately, we defined the relative noise level, i.e., the absolute noise divided by the loop area. The relative noise level for the simulations in Fig. 6 and 7 (absolute range: 0.1–4 nV) ranges from 1.3 × 10−3 to 50.9 × 10−3 nV/m2. The corresponding range in absolute noise levels for a 5-m loop is 0.05 to 2.0 nV and for a 20-m loop 0.4 to 16.0 nV. We chose the noise level of 0.125 nV for the 5-m, 0.5 nV for the 10-m, and 2 nV for the 20-m loop. As above, we repeated the runs for each combination 100 times. The q range was varied depending on ztable to achieve optimum resolution properties at the depth of interest. In doing so, the maximum q to achieve an adequate resolution for the maximum ztable (10 m) was 10 A s for the 5-m loop, 6 A s for the 10-m loop, and 3.2 A s for the 20-m loop. Assuming a maximum pulse length of 10 ms, a maximum transmit power of 1000 A (5-m loop), 600 A (10-m loop), and 320 A (20-m loop) would be necessary to generate the q range of these simulations in reality. Existing devices have a maximum power of 600 to 800 A (Walsh et al., 2014). It should be noted that, at least for a loop diameter of 5 m, an unsaturated zone with 10-m thickness cannot completely be resolved by MRS with currently available equipment.

In Fig. 8 (Model 1) and Fig. 9 (Model 2), we show again median, quartile, outliers, and the square root of the predicted variance to visualize the performance of the WR parameter estimations. At least for Model 1 (Fig. 8), we observe again clear relationships. Except for θS, the accuracy of the estimation decreases with increasing ztable and with increasing loop size. The accuracy of θS slightly increases with loop size and decreases with ztable. However, the accuracy of the θS estimation is much better than the other two parameters, as for the cases discussed above. The uncertainty in θS is <0.025 for all realizations. This is true for both Models 1 and 2, whereas the performance of the other two parameters is very different in both models. For h0 in the case of Model 1, we again observe the trend to underestimation with increasing ztable. The quartile is very similar for each loop size, whereas the risk of outliers increases with increasing loop size. As in the analysis above, the outliers tend to overestimate h0 rather than underestimating it. Consequently, when considering σpred as an uncertainty measure, one must pay attention to the fact that it describes mainly these overestimations. This is particularly true for ztable < 6 m. For deeper regions, σpred overestimates the uncertainties for h0. For the n estimation in Model 1, σpred appears to be larger than the actual uncertainty. This is obvious for the 5- and 10-m loops at ztable > 4 m, while for the 20-m loop, σpred exceeds the plotting range. The quartile for n as a function of ztable is similar for all loop sizes and shows a constant behavior from ztable = 1 to 8 m. The risk of outliers generally increases with loop size and ztable, and again, overestimations of n are more likely than underestimations. For ztable = 1 m, this risk is higher than for ztable > 2 m. Generally, the resolution of MRS in the direct vicinity of the measurement loop is decreased because the sensitivity curves κ(q,z) exhibit local maxima that increase for very small q (see Fig. 2b). Consequently, we observe the decrease in the outlier range for shallow ztable also for Model 2 (Fig. 9c, 9f, and 9i). In contrast to Model 1, we do not observe an increase of the outlier range for increasing ztable for Model 2. Regarding the h0 estimates for Model 2, the quartile is generally larger than for Model 1, and both quartile and outlier range increase with loop size. However, the increase in the outlier range with ztable is smaller for Model 2 than Model 1.


As the synthetic study shows, the performance of the CF inversion does not depend only on several measurement parameters but also on subsurface properties. In general, the accuracy of the WR parameter estimation increases with the number of applied pulse moments and, of course, decreases with noise level. The parameter h0 shows a trend to be underestimated. This tendency rises with increasing noise level and increasing ztable. Furthermore, the accuracy of the PDI estimates increases with increasing width of the CF. For a large CF, the accuracy of the PDI estimates decreases with decreasing ztable, whereas for narrow CF it decreases with increasing ztable. The use of larger loops (i.e., >10-m diameter) cannot improve the PDI estimates. Therefore, we suggest using a measurement loop as small as possible for a reliable CF inversion. However, a small loop has limited penetration depth compared with a larger loop at the same q. Hence, the application of the CF inversion is limited to a maximum depth of about 10 m, while an accurate WR parameter estimation at increasing depth can be expected only for large CFs.

Real Data Example

We next applied the CF inversion to the data set shown in Fig. 1d. Figure 10 shows the results together with the covariance matrices and the uncertainty estimates. The covariance matrix is very similar to those shown in Fig. 4 and 5. The uncertainty intervals indicate that the WR parameter estimations are reliable. In Table 3, the parameter estimates for the three WR models are shown and compared with reference values measured on samples in the laboratory. Four samples were taken in the middle of the measurement loop at about the 50- to 60-cm depth. We assumed that the material from this depth range represents the properties of the entire sand layer, which was found to be lithologically homogenous down to a depth of about 7 m (Costabel and Yaramanci, 2011a). The WR reference measurements for the samples were conducted using the evaporation method (Peters and Durner, 2008) and approximated with the BC, VG, and KO models. The corresponding values in Table 3 represent the parameter estimation ranges of the four repetitions, including their 95% uncertainties. All WR values are in the expected ranges, except for the θS values, which seem to be smaller than usually found for sand (e.g., Hammecker et al., 2004). This is because the samples were saturated in a water bath, and thus a certain amount of residual air in the pore space could not be avoided.

In addition to the MRS from spring 2013 shown in Fig. 10, three further MRS measurements were conducted at the same position, one on the same day with a larger loop of 11-m diameter. The other two soundings were conducted in autumn 2012. Except for the MRS measurement in 2013 with the larger loop, the estimates h0 and the PDI for each WR model are within a reliable range. The θS estimates of the MRSs differ from 0.21 to 0.24 and are underestimated compared with the reference values ranging from 0.25 to 0.29. As expected beforehand, a certain portion of the pore water in the smallest pores was not visible with MRS due to instrumental dead time. This undetectable portion was expected to be similar to θR but it was actually larger. We note that, at least for the example of this study, the MRS water content underestimated the free (i.e., moveable) water content by 3 to 5%.

The WR parameter estimation using the MRS in spring 2013 with the 11-m loop failed. Although θS could be estimated with high accuracy, the estimates for h0 and the PDI had uncertainties of >100%. The reason for the failure is probably the shallow groundwater table in this season (ztable at 1.3 m). As shown with the synthetic data, the shallower the ztable is, the smaller the measurement loop should be. In contrast to the measurement in 2013, the MRS with the 11-m loop from autumn 2012 provided reliable estimates because ztable was deeper (1.8 m). The two MRS measurements with the small loop size of 6-m diameter provided similar results, although ztable was at different depths. This is because the smaller loop resolved the shape of the CF better, which was also expected from the synthetic study.


We investigated a dedicated inversion method for MRS measurements with a focus on the vadose zone, particularly on the CF, which we refer to as CF inversion. The CF inversion directly provides WR parameters under the assumption that the CF is in equilibrium and the substrate is homogeneous. Synthetic and field data demonstrated that θS and both the shifting parameter h0 and the PDI of the WR models of Brooks and Corey (1964), van Genuchten (1980), and Kosugi (1996) can be determined if the water table position is known and included in the inversion as prior information. The CF inversion can provide reliable h0 and PDI estimates with parameter uncertainties slightly higher than from laboratory measurements. Our synthetic study suggests that the CF inversion can be applied if the ztable is located at a depth ranging from 1 up to 10 m. To increase the quality of the parameter estimates under unfavorable noise conditions, the number of stacks or the number of pulse moments can be increased. The latter is often less time consuming in practice. Apart from an adequate number of pulse moments, the CF inversion requires an MRS measurement layout, i.e., particular loop size and q range, that is adjusted with regard to the expected depth and size of the CF.

Our parameter study did not include any systematic uncertainties that might result from possible violation of the equilibrium or homogeneity assumptions throughout the relevant depth range. In practice, discontinuities might occur from seepage water or from lithologic changes inside the unsaturated zone. Consequently, it is not recommended that this inversion scheme will be applied to data measured directly after rainfall or at locations with significant soil inhomogeneity. For instance, the real data set presented here showed a relevant increase in θMRS in the topsoil at z < 0.3 m. However, for this case we showed that after excluding the MRS signals carrying the topsoil information, the CF inversion could reliably be applied. In the future, more experience in applying the CF inversion in MRS is desired to assess its practical application and its sensitivity to the underlying, relatively strong, assumptions. Nevertheless, the suggested approach represents the first attempt to estimate WR parameters noninvasively and inclusively at a plot scale of about 25 to 100 m2, which is the footprint of a small MRS measurement loop.

In general, MRS equipment optimized for vadose zone investigations must provide the ability to generate a high energizing current on demand and, at the same time, the ability to generate short pulses (Walsh et al., 2012, 2014). In this regard, we expect further developments of MRS equipment in the future improving the measurement performance for vadose zone investigations, e.g., the determination of θR. Future research on MRS in this field should focus on time-lapse measurements combined with natural or induced water infiltration. We expect that monitoring the infiltration process with MRS can estimate hydraulic parameters with higher accuracy, at least for vertical water flux. In addition, the whole NMR time series should be taken into account for future interpretation approaches to take advantage of the additional pore information that is encoded in the relaxation times.

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Open Access