Describing hydraulic properties in the subsurface in at least two dimensions is one of the main objectives in hydrogeophysics. However, due to the limited resolution and ambiguity of the individual methods, those images are often blurry. We have developed a methodology to combine two measuring methods, magnetic resonance tomography (MRT) and electrical resistivity tomography (ERT). To this end, we extend a structurally coupled cooperative inversion scheme to three parameters. It results in clearer images of the three main parameters: water content, relaxation time, and electrical resistivity; thus, there is a less ambiguous hydrogeophysical interpretation. Synthetic models demonstrate its effectiveness and show how the parameters of the coupling equation affect the images and how they can be chosen. Furthermore, we examine the influence of resistivity structures on the MRT kernel function. We apply the method to a roll-along MRT data set and a detailed ERT profile. As a final result, a hydraulic conductivity image is produced. Known ground-penetrating radar reflectors act as the ground truth and demonstrate that the obtained images are improved by the structural coupling.

Informing hydrologic modeling with subsurface information obtained from nondestructive geophysical imaging has become an established and broad field of research called hydrogeophysics (e.g., Rubin and Hubbard, 2006). Many geophysical imaging techniques are used to provide structural information on the extent and position of hydrogeologically relevant units as well to investigate their physical parameters to describe the hydraulic properties (e.g., Sulzbacher et al., 2012; Binley et al., 2015; Flinchum et al., 2018).

The most commonly used methods in hydrogeophysics are electrical resistivity tomography (ERT), ground-penetrating radar (GPR), seismic ray tomography (SRT), and nuclear magnetic resonance (NMR). To provide a more comprehensive set of information and to decrease the level of model uncertainty in inversion and interpretation, multiple techniques can be jointly applied to complement each other (e.g., Linde et al., 2006).

There are numerous examples of joint application and joint or constrained inversion of data from ERT and seismic methods (e.g., using seismic refraction data) (Gallardo and Meju, 2003) or GPR (Doetsch et al., 2012), GPR with SRT (Linde et al., 2006), or even the three of them (Doetsch et al., 2010). However, to date the involvement of the emerging surface-NMR method into joint inversion is limited. We focus on coupling the surface-NMR technique and ERT. Contrary to all other geophysical techniques, surface NMR measures the response of the hydrogen proton and enables not only detecting groundwater but also quantifying the amount and distribution of subsurface water by a direct proportionality to the signal amplitude. In addition, the decaying signal provides access to the hydraulic conductivity via proxy relations to the relaxation times T2*,T2, and T1. Several case studies in the past few decades have demonstrated the use of surface NMR for aquifer characterization (e.g., Costabel et al., 2017; Flinchum et al., 2019).

Most surface-NMR applications deal with 1D measurements referred to as magnetic-resonance sounding, i.e., estimating water content and relaxation time depth profiles. Two-dimensional tomography referred to as magnetic resonance tomography (MRT) has been first presented by Hertrich et al. (2007) and allows for subsurface imaging. Recently, efficient layout schemes were introduced (Jiang et al., 2015), using an elongated transmitter loop and an in-loop receiver array (ETRA) promising an increase in the use of MRT. Dlugosch et al. (2014) present the first QT-type inversion (Mueller-Petke and Yaramanci, 2010) using the full decays for MRT measurements enabling access to image the relaxation time distribution. However, this important link to image a 2D distribution of hydraulic conductivity has not yet been applied.

Independent of the dimensionality of the inverse problem, accurate knowledge about the subsurface resistivity distribution is required to calculate the surface-NMR kernel function, especially when considering the complex-valued kernels (Weichman et al., 2000) in a complex inversion scheme (Braun et al., 2005). Various strategies handling resistivity distribution have been presented. Braun et al. (2009) show that resistivity can be estimated from surface-NMR 1D inversion with block discretization. However, this development has not yet been followed further due to the strong nonuniqueness of the solution as shown by Lehmann-Horn et al. (2012). Instead, joint inversion using the transient electromagnetic method (Behroozmand et al., 2012) or vertical electrical soundings (VES) (Günther and Müller-Petke, 2012) has shown practical applicability and, equally important, provides improved resolution and hydrogeologic characterization. Alternatively, Günther and Rücker (2006) study using arbitrary smooth resistivity distributions based on the iterative reweighting of model gradients. Later, using structural constraints and studying coupling parameters in detail, Skibbe et al. (2018) demonstrate that layered and smooth models can be derived from joint inversion of surface NMR and VES and high-resolution models can be obtained.

Braun and Yaramanci (2011) study the impact of 2D resistivity distributions on the kernel function and the NMR sounding curves. Calculating accurate magnetic fields on the kernel discretization based on 2D resistivity structures is, however, not straightforward. Recently, Skibbe et al. (2020) implement arbitrary smooth resistivity distribution (e.g., obtained from ERT inversion) and publish the open-source code COMET based on the modeling and inversion framework pyGIMLi (Rücker et al., 2017). MRT is sensitive to resistivity estimates for amplitude inversion (Dlugosch et al., 2014) and for complex inversion (Skibbe et al., 2020). Complex MRT inversion based on the concepts for 1D complex inversion is presented by Müller-Petke et al. (2016). However, all earlier implementations use simplified geometries with piecewise resistivity (Braun and Yaramanci, 2011; Lehmann-Horn et al., 2011; Jiang et al., 2018, 2020). In this context, the work of Jiang et al. (2020) is interesting because MRT inversion was structurally coupled with GPR and improved MRT images were presented. Until now, no joint inversion of MRT and resistivity has been presented, but, according to the 1D joint inversion results, it can be hypothesized that improved hydrogeologic characterization can be expected.

Whereas for 1D assumptions a joint inversion can be easily obtained by common layer interfaces (Behroozmand et al., 2012; Günther and Müller-Petke, 2012), there are a few different concepts for joint inversion of 2D parameter distributions. Typically, a cooperative inversion aims at common gradients in a prediscretized model discretization going back to ideas presented by Zhang and Morgan (1996) or Haber and Oldenburg (1997). There is the widely applied method of cross gradients introduced by Gallardo and Meju (2003) and applied in hydrogeophysical imaging by Linde et al. (2006). The technique requires computation of the local model gradients and was originally restricted to regular grids. For irregular meshes, Lelièvre and Farquharson (2013) introduce gradient operators, but they point out null-space problems in gradient computation for triangular meshes so that parts of the model can appear decoupled, e.g., for half-regular triangulations. Recently, Jordi et al. (2020) solve the problem of local gradient computation by using geostatistical operators as larger footprints over which the gradients are computed more with improved stability and independent of the mesh discretization. An alternative approach has been presented for 2D ERT and SRT (Hellman et al., 2017; Ronczka et al., 2017) using a weighting of the model boundaries according to the gradient in the other parameter. This method was first extended to the three MRT/ERT parameters using a 1D parameterization.

In this paper, we extend the structurally coupled inversion method to the first application on 2D MRT and ERT as typical hydrogeophysical methods. The mutual benefits of all parameters and the choice of coupling parameters are investigated by simple synthetic models. We study the influence of ERT-derived resistivity on the kernel and analyze the impact on the final inversion results for different schemes of updating resistivity and kernel function during the inversion. Finally, we show a field case in which for the first time the 2D hydraulic conductivity distribution is derived from the coupled MRT and ERT inversion.

Theory

This section provides a short introduction on the applied structurally coupled cooperative inversion (SCCI) method and demonstrates it by using a simple case.

In general, the SCCI algorithm acts, independently on the used geophysical method, by locally lowering the smoothness constraints in the minimization problem. The SCCI works directly with the coupled parameters and calculates a combined smoothness constraint based on the roughness of the coupled parameter distributions, in our case, resistivity, water content, and relaxation times. This means that we use two (or any number of) different methods and can couple the results via SCCI without changing the original inversion algorithms. A prerequisite for the SCCI is that all parameters are defined on the same geometry using the same shaped constraint matrix C0. A comprehensive description of the SCCI coupling parameter in one dimension can be found in Skibbe et al. (2018). Applications in two dimensions, although using other geophysical methods, are given by Hellman et al. (2017) and Ronczka et al. (2017).

There are three main steps in calculating the constraint weights; we outline them here shortly and visualize them in the “Synthetic study” section using a demonstration case. First, after each iteration of the single inversion instances, the roughness vector rp of each parameter distribution p is calculated, resulting in one roughness value per inner mesh boundary that divides two adjacent cells:
rp=C0mp=C0tp(p).
(1)
We prefer to use the transformed (linearized) model mp=tp(p) instead of the physical parameter values p (resistivity, water content, or relaxation times) for equation 1, and the resulting vector r is referred to as normalized roughness. As the coupling exchanges information between the parameters, the transformation also ensures that the roughnesses are in comparable orders of magnitude. For the water content, we use a cotangent transformation within lower and upper bounds of 0%–70%, respectively, as explained by Dlugosch et al. (2014). The resistivity and relaxation time distributions are transformed using a double-logarithmic transformation following Guenther (2004) and Kim and Kim (2011), with boundaries of 0.005–1 s for the relaxation times. Those boundaries are wide limits of present and resolvable parameters in nature. Using the normalized roughness, for each of the n coupled parameters p1pn, an array of individual weights w1wn is calculated:
wp=a|rp|+a+b,for  pin[1n].
(2)
The general behavior of the SCCI equation (equation 2) is shown in Hellman et al. (2017, Figure 2). Note that Hellman et al. (2017) use three SCCI parameters (a, b, and c), but due to earlier investigations (Skibbe et al., 2018) the third parameter c is, however, not necessary and is set to c=1, which then results in our equation 2. The values for the coupling parameters a and b are confined between 0a/b1. The individual weights w1wn are then combined to wc1wcn and multiplied with the constraint matrix C0 of the single inversion instances before the next iteration. As the constraint weight of one parameter, we use the combination of each other parameter’s individual weights (e.g., for the water content constraint weight, we use the individual weights of the resistivity and the relaxation time):
wcp=i=1ipwi,for  pin[1n],
(3)
and thus
Cp=diag(wcp)·C0.
(4)

In addition to the coupling parameters, a global minimum of the constraint weight demonstrably leads to increased stability when coupling more than two parameters. Two effects take place during the SCCI: first, the exchange of structural information between the parameters and, second, the focusing on common boundaries. If two parameters are coupled, then naturally the minimum weight is equal to b. If more parameters are coupled (e.g., three when coupling MRT and ERT), the minimum of the combined weights can be smaller than b, due to the multiplication of weights (equation 3). In those cases, the minimum weight has to be set explicitly; the default value for the minimum of the combined weight is b. This ensures that the focusing effect takes place over several iterations to give the first effect of exchanging information between the methods more time to find the common boundaries of all parameters and not just a blocky representation of each parameter by itself.

Hydraulic conductivity from NMR

The target of a hydrogeophysical investigation can be such a zonation into distinct hydrogeologic units, e.g., represented by cluster center values. However, typically we are more interested in hydraulic parameters than geophysical ones. The water content is already a target parameter of the MRT as is the porosity to a certain degree, when assuming fully saturated media. Taking the relaxation time into account, we are able to compute the most important quantity for hydraulic imaging, i.e., the hydraulic conductivity. There are some simple formulas from borehole logging and more sophisticated ones such as the Kozeny-Godefroy model (KGM) presented by Dlugosch et al. (2013):
KKGM=δg8τ2ηΦ(Dρ+(Dρ)2+4DTBT2TBT2)2.
(5)

The KGM estimates the hydraulic conductivity KKGM based on the surface-NMR parameters’ water content (used as porosity Φ, assuming fully saturated media), the relaxation time T2, the temperature-dependent physical properties’ bulk relaxation time TB, the dynamic viscosity η, the self-diffusion coefficient D, the density ρ, tortuosity τ, and a measure of surface relaxivity δ describing the reactivity of the inner surface. The KGM works for unconsolidated materials in the absence of clay and for laminar flow conditions.

Demonstration case

The coupling usually starts after an initial smoothness-constrained inversion of the single methods is performed. In this case, we use the simple geometry of a circle (Figure 1a) and perform ERT and MRT inversions for a simulated synthetic data set based on the model in Table 1. The layout of the MRT and ERT will be discussed in more detail along with the more comprehensive synthetic calculations later on. We will not discuss the smooth single inversion results of ERT and MRT and instead focus on the 2D visualization and explanation of the SCCI scheme that follows after.

For the smooth inversion images, a constraint weight is calculated using the mentioned equations and visualized in Figure 2a2c. All of the parameters illustrated in Figure 2 share the same color bars (one for each parameter), with the synthetic model marked with white lines for reference. The individual weights wres, wwc, and wt2 are represented as black lines, with the thicker lines representing stronger parameter gradients between the adjacent cells. The gradients of the smooth initial parameter distribution are naturally distributed over a large area and many cell boundaries. At the first iteration, the coupling interacts with nearly every gradient, but the new weights for the constraints will also be smoothly distributed. Only locally, where multiple parameters show gradients, will the combined weights be low enough to introduce changes from one iteration to the next. After the first coupled SCCI iteration, the local lowering of the smoothness constraint eventually leads to a sharpening of the model gradients (and new weights). The combination of all other weights is used to lower the constraint for each parameter (equation 3) to ensure that the gradients are converging to common boundaries for all parameters. The reevaluated results and weights after the first iteration are shown in Figure 2d2f. Usually, several iterations are needed for the models and therefore for the constraint weights to reach a steady state at which the SCCI is finished, indicated by a negligible change in the model and the weights. The final (the eighth iteration) results and the corresponding constraint weights are plotted in Figure 2g2i, leaving only one major parameter contrast that coincides with the synthetic geometry (the thin white circle). Note that the data misfit is comparably equal for all iterations and does not improve during SCCI (Table 2).

The resulting models are more blocky by nature compared with the smoothness inversion, which is expected to be helpful if common parameter gradients of the used methods are to be resolved. Compared with established methods such as cluster analysis, SCCI searches the parameter gradients driven by the data. This leads, compared with other methods, not only to sharper boundaries in two dimensions but also fits the actual parameter values of the different regions of the models better (compare Figure 2 with the values marked in the color bar) because the parameter distributions are changing during the SCCI to find common boundaries instead of mapping different regions.

We split this investigation into two parts using the same measurement setup. We already introduced a simple target structure in Figure 1. The first part of the synthetic study will focus on the use of the SCCI parameters a and b. Therefore, we deliberately chose a simple nongeologic model with strong NMR parameter contrasts and a resistivity value that we know has an effect on the NMR kernel function. This serves the purpose of validating the algorithm and showcasing the effects of the coupling as well as the resistivity later on.

In preparation of the structural coupling, separate standard smoothness-constrained inversions are performed for ERT and MRT. To cover both presented target structures with the ERT measurements, we use 121 electrodes in a 1 m spacing from −60 to 60 m. A multiple-gradient array is used for the ERT because it provides a good trade-off among accuracy, depth-of-investigation, and speed because it is especially designed for the use of multichannel devices. We model a total of 2842 synthetic ERT data points and add Gaussian distributed noise with a standard deviation of 2% of the data. For MRT, we chose a layout of three 40×40m quadratic loops, which are positioned, half-overlapping, between −40 and 40 m. We use all possible combinations of the loops leading to nine different configurations and use 24 logarithmically distributed pulse moments ranging from 0.28 to 13.56 As. The earth magnetic field is assumed to have a declination of 90° and an inclination of 60° with the loops’ profile direction pointing toward the positive x-axis, the depth is defined negative downward, and the magnitude of 48,000 nT leads to a Larmor frequency of 2043 Hz. A Gaussian-distributed random noise with a standard deviation of 10 nV is added to the synthetic data set (real and imaginary separately) before the gating of the data.

The smooth inversion results of the ERT and MRT used as starting models for all following studies are shown in Figure 3a3c for the three parameters of resistivity, water content, and relaxation time, respectively. For reference, the synthetic model values are marked in the color bar for reference and the shape of the anomaly is displayed in thin white lines. Note that these lines are not included in the meshes.

SCCI is performed for several iterations until reaching a steady state. Overall, we compute eight iterations for both methods to reach a steady state. Note that all iterations retrieve results with equivalent data misfit as already obtained by the smooth inversion. On the contrary, notable changes in the data misfit are only observed if either the chosen SCCI parameters are too small and/or if the smooth inversion was already over- or underfitted.

SCCI parameters

This part of the synthetic study aims to determine the effect of the coupling parameters on the 2D images. We investigated a wider and finer range of values for the parameters a and b, but we will focus the discussion on three values (nine combinations). The values for a and b include 0.01, which already reacts on very small gradients and potentially leads to overfitting and other problems. The second value of 0.1 has proven to be a reliable choice (Hellman et al., 2017; Skibbe et al., 2018), whereas the third value of 0.3 only leads to small changes of the constraint weights as we will show next.

The final SCCI result for the nine combinations is shown in Figure 4 for the relaxation time distribution as an example. The corresponding images for resistivity and water content times can be found in Appendix  A for completeness. In general, they behave similar to the relaxation time.

From the used equation 2, as well as the experiences of other studies involving SCCI, several observations can be made. The choice of parameter a determines the threshold at which a parameter contrast triggers a lowering in the constraint weight. Low values of a will allow even the smallest contrast to lower the constraint weight, whereas large values will allow the strongest gradient to have this effect. Inspecting the first row of Figure 4 shows that the smooth nature of the starting model is kept throughout the SCCI and results in little change of the general model behavior. This effect is independent of parameter b for small values of a.

Parameter b is adding a general value to the constraint weight, thus stabilizing the constraints. In addition, the minimum of the combined constraint weight is set to b by default. For small values of b (column 1), the parameter a has total control over the final result, whereas for large values of b (column 3), the effect of a is practically nonexistent.

The influence on the coupling parameters becomes even more clear when looking at the constraint weights directly instead of on the parameter distributions. Figure 5 shows the combined constraint weights wct2 following equation 3 for the smooth relaxation times distribution of Figure 3c and different coupling parameters. As for the demonstration case, thicker black lines indicate that according to equation 3 the parameter contrast leads to a lowering of the smoothness constraint for the other two parameter distributions in the next iteration of SCCI. Areas without black lines will not be influenced directly by SCCI, meaning that the original smoothness constraint is not altered. That parameters a and b have an influence on the final SCCI result as can be seen in the different first constraint weight distributions.

The final constraint weights are shown in Figure 6, now calculated for the nine distributions shown in Figure 4. Depending on the initial weights, different features of the starting model are enhanced or neglected during SCCI. However, it is interesting to mention that the synthetic geometry of the layer can be identified for most of the chosen values of a and b (if a is not too small).

Similar to Skibbe et al. (2018), we compare the final SCCI models with the synthetic models by calculating the root mean square (rms) of the models in the displayed inversion area. This gives us a measure of how the different (physical) parameters respond to SCCI. The results for the full range of investigated coupling parameters are shown in Figure 7. In general, larger rms values can be observed if both parameters are small, or if b is too high. This is expected behavior: b is stabilizing the weights; if the value is too high, all gradients are flattened and no common boundary is found. On the contrary, if both values are too low, every gradient leads to a lowering of constraints and, again, no common boundary is established.

The lowest/best rms values are achieved for values of approximately 0.1 for both parameters. The three parameters, or two methods for this regard, behave slightly differently during SCCI. This is probably due to the model transformations of the parameters, and the enhanced gradients are treated differently in the single inversions.

Influence of resistivity

The first part of the synthetic studies hold some general statements concerning SCCI. Apart from the SCCI parameters, other coupling effects may occur; in the special case of ERT and MRT, the resistivity distribution is not only coupled via the SCCI algorithm but is also used for the sensitivity calculation of the MRT. In this section, we discuss to what degree the resistivity can and should be incorporated in the inversion process when coupling ERT and MRT.

In all of the discussed cases, the ERT inverted resistivity distribution has been used in the kernel calculation of the MRT smoothness inversion. In the first case, the SCCI is performed coupling only the MRT parameters’ water content and relaxation times; thus, the resistivity only has an influence on the kernel function but not on the coupled inversion results. The final results using a suitable set of SCCI parameters are displayed in Figure 8a8c.

Further incorporating the resistivity, the SCCI is also performed coupling the resistivity with the other two parameters in the second case. This not only changes the resistivity, but it also influences the final distributions of the water content and relaxation times (Figure 8d8f) through the coupling. All three parameter distributions are better resolved compared to when they are uncoupled from the ERT.

The most comprehensive use of the resistivity would be to recalculate the kernel function for MRT after each iteration of SCCI (the third case). However, this is a time-demanding thing to do. The final results of the third case are shown in Figure 8g8i but show little to no difference compared to the second case (Figure 8d8f).

For reference purposes, we also performed SCCI with kernels based on the synthetic (real) resistivity; this way, the results are only influenced by the coupling and not through the initial kernel function.

The used kernel functions for the coincident case of the first loop are plotted in Figure 9, with respect to the used resistivity distributions, including the smooth start model, the eight upgrades of the fully coupled SCCI, and the synthetic (real) resistivity (Figure 9a9c, respectively).

Although the shape and values of the anomalous region in the half-space can be resolved by SCCI, it is interesting to notice that the behavior is different when the resistivity is changed (compare Figure 8a8c with Figure 8d8f). The smooth resistivity seems to resolve the anomaly in the subsurface sufficiently for the purpose of a kernel calculation because both kernel functions (with or without resistivity update during the inversion) are close to the kernel based on the synthetic resistivity, which means that there is not much room for further improvement. As expected, even with continuous kernel upgrades, no further benefit can be observed in the SCCI results (compare the second and third cases).

To show the performance of SCCI under field conditions, data from the well-known Schillerslage test site near Hannover, Germany, are used. Following Jiang et al. (2020), we characterize the subsurface geometry of the site by an alternating succession of till and meltwater sand and gravel, overlain by aeolian deposits. The meltwater deposits act as aquifers, and the till beds act as aquitards.

Previous, multiple geophysical surveys, boreholes, and laboratory analyses provide the opportunity to verify our results. The MRT data first have been used by Jiang et al. (2018), but only a small subset (one out of three roll-along layouts). Jiang et al. (2020) present a more comprehensive study using the complete data set including GPR, ERT, and borehole measurements, although with the objective of using GPR reflectors for constraining the MRT image.

We use the same MRT and ERT data as described in Jiang et al. (2018, 2020) to conduct our SCCI, but we use the GPR boundaries as a reference only. In the subsurface, we have a shallow bounded aquifer down to a depth of approximately 20 m, and it is separated into two parts by a thin till layer. In the imaging, we focus on the first shallow aquifer only because the second is at the resolution limit and can only be detected by extreme anisotropic regularization and/or meshes adapted to the thin bed. The groundwater table is at a depth of approximately 2 m, and it is expected to be the layer boundary with the strongest contrast in all parameters. The second boundary is the base of the first aquifer, which is undulating between 6 and 14 m along the profile. The aquifer itself consists of several periglacial deposits from fine sand over medium sand to gravel, which is present in pockets of the till surface.

The measurement setup consists of three elongated in-loop receiver array layouts (Jiang et al., 2015) as roll along, resulting in a total of 22 half-overlapping 20×20 receiver loops and three edge-to-edge positioned 60×20 transmitter loops between −55 and 125 m. For all measurement configurations, we used between 20 and 24 logarithmically distributed pulse moments ranging from 0.05 to 9.5 As with an average of 50–60 stacks.

In addition, we use the data from an ERT setup using 240 electrodes between −105 and −173 m with an electrode spacing of 1 m in the middle part and 2 m spacing for the 20 electrodes at each side of the profile. A total of 16,932 data points were recorded using a multiple-gradient array.

As a prerequisite for the MRT inversion, we performed a smoothness constraint inversion of the ERT. Following the conclusions from the investigation concerning the influence of the resistivity on the MRT kernel and inversion, we decide to use the smooth ERT result for the magnetic field calculation for the kernel functions of the MRT. Figure 10 shows the derived smoothness constraint MRT and ERT inversion results. The dotted white lines represent the GPR reflector boundaries found in the investigation of Jiang et al. (2020) and are for validation purposes only; they have not been used in any step of the mesh generation or inversion processes.

Figure 11 shows final inversion results of the SCCI again together with the GPR reflectors in white as a reference. We focused the view on the common resolution area of both methods; however, the inversion mesh is large enough to include all electrodes and the full kernel function of the MRT to eliminate any artifacts from the mesh boundaries.

The vadose zone is characterized by resistivities above 1000Ωm, a water content of 10%, and relaxation times of 100–150 ms. The main contrast of the groundwater table to the first aquifer is clearly resolved in all distributions and coincides with the findings of previous investigations. Due to additional borehole information, it is known that the aquifer consists of gravel marl, supposedly with varying degree of sorting. Higher values for water content (up to 35%) and relaxation times (up to 150 ms) at the profile meters −50 to −25, 0–50, and 75 m onward suggest a better sorted gravel than in the other areas. This actually agrees with the findings of Jiang et al. (2020), although the high vertical constraints they are using might be smoothing this effect, where the SCCI seems to resolve this more clearly. The resistivity distribution shows fewer lateral changes with values of approximately 300Ωm. The lower boundary of the first aquifer follows the GPR reflector in general with deviations at the end of the profile, due to resolution problems. Below the aquifer, an aquitard with low resistivities of 30Ωm and low water contents of approximately 10% is observed, which fits the expected till found in boreholes. The relaxation times show lateral variances between 80 and 200 ms.

Generation of hydraulic models

The target parameters of a hydrogeophysical investigation are often not the geophysical parameters but hydraulic properties that can be derived from them. With the available information, we are able to generate a hydraulic model using the KGM shown in equation 5. Due to past investigations at the field site, we have reliable estimates for the surface relaxivity of 130μms1 and standard values for the tortuosity of 1.5. Values for the physical properties of water, TB, η, D, as well as ρ are derived using a middle temperature of 10°C similar to what is presented in Dlugosch et al. (2013). The KGM is defined for values of T2 and not T2*, however, in this specific case, borehole data confirm that these relaxation times show an agreement for the inversion area (Jiang et al., 2020, Figure 9). For this case, it is possible to substitute T2 with T2*.

Figure 12 shows the hydraulic conductivity distribution based on the final SCCI result of the inverted field case and the aforementioned estimates for the physical properties of the subsurface. The lateral variation of the aquifer system is also seen in the hydraulic conductivity. This fits with the findings of borehole measurements and sediment logs presented in Jiang et al. (2020) showing lateral heterogeneities along the profile. Because we assume better sorted, coarse material where the water content was greater, higher values of hydraulic conductivity make sense. In this case, we observe values between 5×104ms1 and 1×103ms1 as peak value, whereas the aquitard drops below 1×105ms1. Note that fully saturated media are a requirement when the KGM is used; therefore, the first meters are not an adequate representation of the real situation. Therefore, we mask the values in the corresponding Figure 12. With the various prerequisites, there are several uncertainties that need to be accounted for to use the hydraulic model.

The introduced method of structural coupling produces images of resistivity, water content, and relaxation time that are more similar in terms of the parameter gradients in the model. In all shown cases, it is easier to discriminate different zones of the subsurface, compared to the classical smooth inversion where typically only smooth transitions are shown. The biggest improvement was always observed for ERT, where the smooth inversion transitions were always too deep and subsequently shifted upward by SCCI. In total, the SCCI results allow an improved abstraction of the subsurface into zones of constant parameters. For the two rather simple synthetic models, the results represent a rough-surface variant of the input models with almost constant values.

Similar to what is presented by Hellman et al. (2017), the three parameter distributions could be combined into one image by using a cluster-based model integration that is expected to be even better with three parameters compared with only two. The clear model boundaries of the SCCI result will ease the algorithm to find contiguous zones that are represented by single values of resistivity, water content, and relaxation time.

For the field case, the subsurface structure is more complicated as a result of the geologic deposition with different grain sizes (for detailed interpretation, see Jiang et al., 2020). Comparison with the GPR reflector of the aquifer base (middle line in Figure 11) clearly shows that it is well matched with the largest gradients of all three parameters, at least in the center of the profile. At the outer parts, there are deviations due to resolution limits of the MRT.

We show the results using a = 0.15 and b = 0.1 for the SCCI parameter. Considering the number of geophysical parameter jumps we expect, we chose this set of parameters because it provides a good idea concerning the different lithologic units, which in turn serves our subsequent hydrogeophysical interpretation.

Uncertainties of the hydraulic model

The generated hydraulic model shows that it is possible to generate hydraulic conductivity images from geophysical methods, but they will be subject of further research because several uncertainties need to be addressed when using KGM based on surface-NMR measurements. First, T2* needs to be close to T2; i.e., the earth’s magnetic field has to be homogeneous at the field and pore scales. This can only be investigated by additional laboratory or borehole measurements. For our case, this can be seen in the mentioned sediment logs presented by Jiang et al. (2020). Second, only through additional knowledge do we know that fine-grained materials of short relaxation times cannot be accurately described by the surface NMR. In those cases, the water content is underestimated, which leads to incorrect porosities, one of the input parameters of the KGM. Third, the KGM model is only valid for fully saturated material.

Beyond the transformation to hydraulic conductivity based on the KGM equation, surface NMR and ERT measurements have their inherent ambiguities, leading to uncertainties of the hydraulic model. These ambiguities, however, can be reduced using additional measurements or other information, not necessarily in two dimensions. Jiang et al. (2020) show that better structured parameter distributions are more accurate when compared with borehole data. Instead of predefining these boundaries (the case if there are no available GPR data), we achieve the structured parameter distributions using the SCCI and improve the primary parameters and subsequently the hydraulic model. Furthermore, additional 1D investigations, if possible involving boreholes, are helpful to calibrate the KGM calculation, which in turn offers a straightforward way to interpolate the target parameters between the available logs in two dimensions, based on the geometry derived from the 2D measurements.

Choice of coupling parameters

In general, the coupling changes the parameter distributions with respect to common boundaries while maintaining the data fit of the initial smooth inversion results. That means that, independent of the coupling parameters a and b, the final result of a conducted SCCI will be a valid representation of the subsurface with respect to the inversion evaluation criteria (with a few exceptions).

To narrow down the choice of parameters, to obtain a model fitting the general expectations of an investigated case, we summarize our findings on the influence of the two coupling parameters in general:

  • If both values are small (≤0.01), the final result behaves like a blocky model in one dimension when the number of layers is large. For extremely small values (<<0.01), this can lead to instabilities and models that have no physical meaning, similar to not choosing a constraint at all (e.g., as shown by Skibbe et al., 2018).

  • If the values for a and b are both large (>0.3), the global smoothness constraint is barely reduced and the final result will be close to the result of the smoothness constraint inversion. Note that, by default, the constraint cannot have values larger than 1; thus, no increase of the smoothness constraint is intended.

  • It is mostly beneficial to choose a small value for b (0.05b0.15) to ensure a certain stability and then search for the best value for the parameter a.

  • There is no identical optimum for all parameters because these are case-dependent; however, keep in mind that the coupling affects the linearized models (equation 1), meaning if the used model transformations for the inverted parameters are able to linearize the inverse problem, a satisfactory value for a can usually be found between 0.05 and 0.3. If multiple and/or small parameter jumps are expected, smaller values for a are needed. Greater values of a are used, if large, isolated parameter jumps are to be resolved.

Influence of resistivity

According to the results presented in Figure 8, only a few changes are found in the kernel function (Figure 9) and the inversion results. Obviously, the smooth resistivity result, in our case, is already close enough to the final result to describe the behavior of the magnetic field sufficiently. MRT is not primarily sensitive to resistivity structures; therefore, the smooth result of an ERT inversion is sufficiently capturing the important features leading to accurate kernel functions for a subsequent MRT inversion. Because a full kernel computation takes a certain amount of numerical effort, and it is apparently not needed in between the SCCI steps, we recommend the following procedure: (1) perform a smooth ERT inversion, (2) compute the kernel based on smooth resistivity, (3) run the MRT smooth inversion, and (4) run the SCCI inversion for a fixed value of b and a couple of values for a (v, optionally), then recompute the kernel during the SCCI at your convenience (e.g., whenever structural changes in the resistivity distribution are observed). Nevertheless, there might be cases in which an update of the kernel function during SCCI is necessary in each iteration; however, it was not necessary for any of the cases that we investigated.

Further improvement

The results from the gradual decrease of constraints using the data-driven SCCI are not as good as those achieved by including constraints from GPR reflectors done by Jiang et al. (2020), where boundaries are ignored in the inversion by setting the constraint weight to a fixed value of zero at certain boundaries. For zonation purposes, it is beneficial to include these additional zero-valued constraints; therefore, GPR should be part of the hydrogeophysical workflow whenever possible. However, the measurement conditions differ between GPR and surface-NMR measurements, and the depth of investigation may be different due to subsurface conductivity structures. Both constraint types can be combined and should result in a subzonation, e.g., a rough aquifer zonation by GPR and differentiation into several deposits by MRT. Borehole data, e.g., NMR logging, grain-size analyses, and direct push measurements, could further improve the images at depth by introducing more fixed constraints in the inversion, adapting the model transformation, or using them as a reference model in the inversion.

The choice of common boundaries by the SCCI is to a large degree dependent on the chosen mesh discretization, i.e., fine enough to describe the subsurface and the kernel structures. The created directions of the mesh boundaries will be important for the resulting images. It also stands for debate how the SCCI constraints can be combined with anisotropic regularization. To become more independent of the chosen mesh, geostatistical constraints as used by Jordi et al. (2018) are the better choice. Those geostatistical operators also can be used to compute and minimize model gradients in a cross-gradient joint inversion scheme (Jordi et al., 2020). By applying this computationally more intense method, we expect improved and more mesh-independent results. Furthermore, clustering approaches could be incorporated to further clarify the images.

We adapted an existing method of SCCI to three 2D distributed parameters: resistivity, water content, and relaxation time. The method overcomes the blurriness inherent in the ambiguity of the inverse problem of the individual ERT and MRT methods. For the synthetic cases, this results in clear images that are close to the synthetic models. Different from our expectations, the resistivity image has a small effect on the MRT kernel as long as the general structures are described. For the field case with far more complex structures, the SCCI inversion also improves the images and leads to an improved aquifer zonation and parameterization. As a result, we produced 2D models of the hydraulic conductivity, one of the main target properties. The combination of MRT with ERT in a structural manner, and if possible with GPR as additional constraints and borehole data, is recommended as a workflow for medium-scale hydrogeophysical investigations. The 2D nature of the results is beneficial for further investigations, especially if compared with hydraulic models derived from boreholes in one dimension as the ground truth.

This research was supported by the Deutsche Forschungsgemeinschaft (German Research Foundation) under grant MU 3318/3-1. We thank the two anonymous reviewers and the editorial board for the work they put into the paper during the review process. We also would like to express our gratitude for the continuous work with C. Jiang, during the original field campaign as well as for the many fruitful discussions afterward.

Data associated with this research are available and can be obtained by contacting the corresponding author.

APPENDIX A

Resistivity and water content distributions of the SCCI parameter study

Figures A-1 and A-2 show the inversion results of the synthetic data (the layer case) for the parameter distributions’ water content and electrical resistivity for different SCCI parameters. The structural characteristics are similar to those of the relaxation time shown in Figure 4: higher values for the coupling parameters lead to smoother images, whereas small values lead to structures beyond the resolution limit of the methods. The similarities of the three distributions are based on the fact that the constraint weights are calculated based on all parameter roughnesses. Both parameter distributions are able to detect the layer in shape as well as absolute values (the original values are marked in the color bars) for most of the used coupling parameter combinations.

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.