Structurally coupled cooperative inversion of magnetic resonance with resistivity soundings

INTRODUCTION

Geophysical measurements have proven to give valuable noninvasive information for the characterization and interpretation of hydrologic environments. One of the methods particularly suitable for this task is magnetic resonance sounding (MRS) because it allows for retrieving information about the porosity and hydraulic conductivity (e.g., Lachassagne et al., 2005) based on the primary parameters’ water content and relaxation time. In MRS, large surface coils are used to transmit electromagnetic excitation pulses stimulating response signals from hydrogen protons that are received by the same coil. To extract the demanded porosity and hydraulic-conductivity information from the received signal, exact knowledge of the electromagnetic field propagation is essential. However, an exact calculation of the magnetic fields depends on the knowledge of the resistivity distribution. For example, the importance of sufficient information on deep resistivity structures, especially in the presence of a deep conductive layer, has been shown by Behroozmand et al. (2012a, 2013). On the other hand, Lehmann-Horn et al. (2012) demonstrate that smaller, but shallower, resistivity anomalies can also have a significant influence on the MRS response. Even though it has been shown that the resistivity may be estimated from MRS itself (Braun and Yaramanci, 2008), it is commonly recommended to retrieve resistivity information from other geophysical measurements. Different frequency-domain electromagnetic (FDEM) and time-domain electromagnetic and direct-current methods are capable of providing resistivity information. The most commonly used methods are vertical electric sounding (VES), transient electromagnetics (TEM), and FDEM, either surface based or airborne. All of these methods have already been used to support the inversion of MRS data, each having their own advantages and disadvantages in terms of resolution, depth of investigation, and their behavior in the presence of conductors or resistors.

More recent developments, however, target not only the use of the resistivity structure for the magnetic-field calculation but also coupling the resistivity and MRS data inversions for mutual support. The presented joint-inversion approaches of MRS with TEM (Behroozmand et al., 2012a) or VES (Günther and Müller-Petke, 2012; Akca et al., 2014) clearly show improved resolution and more stable results in the resistivity, water content, and relaxation time distributions. Moreover, ambiguities compared with single inversions are reduced (Günther and Müller-Petke, 2012). To date, all joint approaches for MRS are based on a block/layered discretization following either a Marquardt-type algorithm (Behroozmand et al., 2012a; Günther and Müller-Petke, 2012) or a global inversion scheme (e.g., Akca et al., 2014). To our knowledge, there is no joint approach for inverting a smooth distribution with resistivity soundings. Unlike block inversion, in which the number of layers is fixed and the thickness of the layers is estimated along with the target parameters, for smooth inversion, the thickness of a larger number of layers is fixed and only the target parameters are inverted. Both inversion types have advantages and disadvantages. On the one hand, the block inversion greatly depends on the starting model; thus, it needs a priori information, such as the number of layers and is mostly restricted to 1D problems. On the other hand, smooth inversion often has problems resolving smaller parameter variations or large parameter contrasts because the inverse problem is ill-posed and needs regularization via a stabilizer function. However, the smooth inversion is easily extendable into 2D/3D problems for which it is the only solution. Because 2D problems are envisaged for future research, coupling smooth inversion of resistivity and MRS is the first step in this direction.

The standard stabilizer function for smooth inversions is the $L_2$-norm; i.e., the quadratic norm of the roughness is minimized. Recently, Grombacher et al. (2017) introduce the minimum gradient support (MGS) as an alternative stabilizer function for MRS and compare MGS with standard smooth $L_1$- and $L_2$-norm inversion. The MGS method usually leads to blocky features in smooth inversions, but without the exchange of information between the smoothness of the two target parameters. We follow the idea of structural coupling, i.e., searching for similar models, as presented in Haber and Oldenburg (1997). They introduce the idea of exchanging information between two parameters that have little or no physical relation; only similar structures in the parameter gradients are assumed. An implementation for a structurally coupled cooperative inversion (SCCI) was outlined by Günther and Rücker (2006). Later on, a more robust function to couple MRS with VES/IP 1D models along with cluster analysis to obtain blockier subsurface layers was presented (Günther et al., 2010). However, this simple coupling study neglected the decay time and assumed the resistivity to be known beforehand, thus not solving the full nonlinear problem. This coupling strategy was refined and applied to 2D electrical resistivity and traveltime tomography by Hellman et al. (2017) and Ronczka et al. (2017). This approach results in blockier models compared with the smooth inversion, taking into account the different sensitivity and roughness distributions of the individual methods. We use this coupling approach for a joint inversion of MRS data together with resistivity soundings to obtain smooth distributions of resistivity, water content, and relaxation time, but also allowing for blocky features. Although the produced features are similar to that of MGS, we use a weighted $L_2$-norm as stabilizer function, similar to an L1 inversion. SCCI combines the advantage of a coupled inversion with the enhancement of blocky features.

For the MRS inversion, we invert the whole data cube as a function of pulse moment $q$ and time $t$, also known as qt inversion (Mueller-Petke and Yaramanci, 2010) for water content and relaxation times simultaneously. We reduce the relaxation-time model space to monoexponential relaxation, as presented by Dlugosch et al. (2014). An overview of the currently used inversions, including specifications and classification, can be found in Müller-Petke et al. (2016, Figure 4).

The paper is organized as follows: First, we introduce the weighting function and give a brief review on all parts of the inversion affected by the structural coupling. Second, the overall (separate) inversion scheme for resistivity, water content, and relaxation times is discussed, using a weighting between these distributions based on the actual roughness as a constraint. Third, a detailed investigation of all parameters is performed by using two synthetic models: a smooth model and a blocky model. Next, field data are used to demonstrate the performance and stability of the presented algorithm using a well-known, but complex, site with sufficient ground truth. Finally, the applicability of the structural coupling is discussed in general and with respect to the choice of the coupling parameters.

METHODOLOGY

We introduce only the fundamental equations of the MRS forward modeling. For the MRS theory, we refer to Weichman et al. (2000) or Hertrich (2008). An overview of the state-of-the-art signal processing is given by Müller-Petke et al. (2016). Unlike nuclear magnetic resonance (NMR) applied in a laboratory framework, MRS is mostly limited to the free-induction decay experiment. Thus, in the following, the term relaxation time exclusively refers to $T_2^{*}$. The coupled inversion targets to combine water content and relaxation time with a resistivity distribution, e.g., from VES, for which only a brief introduction is given. However, any electrical or electromagnetic sounding method can be used likewise. The second focus of this section is to introduce the structural coupling inversion approach.

MRS forward calculation

MRS uses the response of the hydrogen proton to an artificially generated electromagnetic pulse. Depending on the pulse moment $q=\tau ·I$ with the pulse duration $\tau$ and a current amplitude $I$ applied to a surface loop, the resulting magnetic field causes an excitation of the protons’ spin orientation. After the pulse, the protons relax to an equilibrium state that is defined by the static magnetic field. The magnetic field caused by this relaxation usually induces a voltage of a few hundred nV in the receiver loop(s).

Assuming a monoexponential relaxation time in each depth layer, the general formulation of the complex MRS signal envelope reads (e.g., Dlugosch et al., 2014)

For each pulse moment $q$, we measure a signal $V$ as a function of time $t$. The signal envelope is characterized by an exponential decay described by the relaxation time $T_2^{*}$, a volume integral over water content $\mathbf{\theta }(\mathbf{r})$, and a sensitivity function $K(\mathbf{r},q)$, referred to as kernel.

The basic formulation of the NMR forward kernel is a function of space $\mathbf{r}$ and pulse moments $q$ for a coincident loop and reads (Weichman et al., 2000)

with the Larmor frequency ${\omega }_L$ and the net magnetization in the equilibrium state $M_0$. The Larmor frequency $f_L={\omega }_L/(2\pi )$ with $w_L=|\gamma B_0|$ depends on the gyromagnetic ratio $\gamma$ and the earth’s magnetic field $B_0$. The sine term describes the spin excitation that depends on $q$, $\gamma$, and the transmitter magnetic field component perpendicular to the earth’s magnetic field and corotating with the spin ${\mathbf{B}}_{+}^{\perp }(\mathbf{r})$. An elliptical decomposition scheme can be found in Hertrich (2008). This sine term causes the strong oscillating behavior of the kernel function. Phase shifts $\xi (\mathbf{r})$ occurring between transmitter/receiver coil and point $\mathbf{r}$ are described by the exponential term. Finally, ${\mathbf{B}}_{-}^{\perp }(\mathbf{r})$ is the receiver magnetic-field component perpendicular to the earth’s magnetic field and counter rotating with the spin. Due to the oscillation of the excitation pattern, a high mesh discretization in the vicinity of the loops is needed to obtain an accurate representation of the kernel function. Because the kernel function must be calculated in three dimensions, this can easily result in meshes with millions of nodes. However, calculating the relatively smooth magnetic field on this kernel mesh directly is neither recommended nor needed to get a sufficient representation of the magnetic field. Therefore, we calculate the magnetic field on a second mesh suitable for interpolation with high accuracy but with fewer nodes than the kernel mesh.

Even though the kernel calculation is always in three dimensions, only 1D inversions are performed. Thus, it is reasonable to confine the resistivity changes to 1D as well. In this case, a semianalytical solution for the magnetic field can be found using Bessel transformations for horizontal electric-dipole (HED) sources. The basic formulations for HED fields and recursion formulas (for 1D resistivity) can be found in Ward and Hohmann (1988). We use Hankel factors according to Christensen (1990) for solving Bessel integrals and provide magnetic fields for the transmitter. To provide sufficiently accurate fields for the kernel function, the dipole discretization can easily reach a couple of hundred dipoles per loop. Even with a semianalytical solution, this becomes a computational bottleneck in the overall calculation of the MRS forward response.

We calculate the first kernel function based on the initial smooth inversion result of the VES data. Further updates of the resistivity distribution are small, and a recalculation of the kernel function at each iteration significantly increases the running time of SCCI, due to the smooth parameter distribution. The influences of resistivity on MRS have been investigated by, e.g., Braun (2007), Braun and Yaramanci (2008), and Behroozmand et al. (2013). To save calculation time, a recalculation of the kernel function can be omitted in cases in which the resistivity update is small or for small loop sizes because the initial kernel is already sufficient for parameter estimation. However, a method to efficiently cache most of the calculation steps and using interpolation to reduce the number of dipole calculations is presented in Appendix A.

Note that Behroozmand et al. (2012b) also present a fast MRS forward calculation. However, they use a few layer block discretization and thus the number of recursions is greatly reduced resulting in a significant speed-up.

Similar to the magnetic field, the VES forward calculation can be solved semianalytically. We use the 801 point filter of Anderson (1989) to solve the Bessel integrals via Hankel transformation. The aforementioned literature concerning Hankel formulations also covers VES forward calculations.

Inversion

The inverse problem can generally be written in the form of a minimization of a data residual

where $N$ is the number of data points $d_i$, $f_i(\mathbf{m})$ are the model responses, and ${\epsilon }_i$ is the associated data errors. In the case of an optimum data fit, the typical evaluation criterion is ${\chi }^2={\mathrm{\Phi }}_d/N=1$. Due to the ambiguity of the inverse problem, an infinite number of models achieves ${\chi }^2=1$. To get a (hydro-)geologically reasonable model, constraints are used to ensure a control of the model characteristics.

The first constraint that we use is a transformation applied to the model vector. For each method, the model vector $\mathbf{m}$ consists of $M$ parameter distributions $\mathbf{p}$, including transformations $t$:

In the case of electromagnetic methods, the model is represented by the resistivity distribution $\rho$, whereas in the case of MRS, the model vector contains the parameter distributions for water content $\theta$ and relaxation times $T_2^{*}$. We use two-sided logarithmic model transformations (Günther, 2004; Kim and Kim, 2011) for resistivity and relaxation times with bounds of 1–10,000 Ωm and 5–1000 ms, respectively, whereas for the water content, a cotangent transformation (Dlugosch et al., 2014) is used, confined between 0% and 70%. Additionally, we append a logarithmic transformation on the data space of the VES.

We solve the nonlinear inverse problem by a Gauss-Newton method calculating model update $\mathrm{\Delta }\mathbf{m}$ in a linearized subproblem. A popular method of inverting 1D parameter distributions is block inversion using a Marquardt-Levenberg scheme. The inversion scheme is used to estimate a fixed number of thicknesses and their parameter. This very small model space, usually just a few layers, leads to an over-determined problem to solve and does not need additional regularization. However, a good starting model is needed to find the global minimum of equation 3 to get a reasonable model for subsequent interpretation.

We use another approach, in which the layer thicknesses are fixed and only the parameter values are inverted. This scheme is often referred to as Occam inversion (Constable et al., 1987), and it usually allows for a lot of parameters, and, therefore, it has a larger model space compared with a Marquardt inversion. This leads to an under-determined inverse problem, which needs to be regularized. The basic minimization for solving the inverse problem (equation 3) now consists of a data $({\mathrm{\Phi }}_d)$ and a model residual $({\mathrm{\Phi }}_m)$ and states as

A regularization parameter $(\lambda )$ between data and model is introduced. The model residual is defined as the $L_2$-norm of model and constraint matrix $\mathbf{C}$:

The constraint matrix is used to ensure a certain behavior of the model. The simplest form is a first-order derivative matrix that is for a 1D discretization:

where the number of columns equals the length of the model vector and the number of rows equals the number of boundaries (the length of the model minus 1).

The value for $\lambda$ has to be chosen in a way that satisfies ${\chi }^2=1$. Additionally, according to Constable et al. (1987), we search for the highest $\lambda$ (the smoothest model) that explains the data within the error bounds. This value is further referred to as optimal lambda ${\lambda }_{\mathrm{opt}}$.

All computations are implemented in the open-source Python framework pyGIMLi (Rücker et al., 2017). We also use the implemented brute-force approach for calculation of the VES Jacobian. The inversion uses an iterative least-squares solver with an inexact line search (Günther et al., 2006) to solve the inverse subproblem. For the MRS inversion, we invert the complex data set, i.e., the real and imaginary parts, independently, as described by Müller-Petke et al. (2016).

Structural coupling

In the presented case, we use two different methods (VES and MRS) inverting for three different parameter distributions (resistivity, water content, and relaxation time) that are coupled with each other. The coupling is applied between the parameters, not the methods. The structural coupling of different parameters is achieved by mutual weighting of the derivative.

First, the roughness of each parameter distribution $\mathbf{p}$ is calculated by

Second, a parameter-specific weighting function $w$ for each boundary $i$ is calculated according to Hellman et al. (2017) by

where $a$ and $b$ are the parameters to control the coupling strength. Additionally, the weights are limited at 1 and are therefore confined between $b$ and 1. Basically, $a$ is used to specify a threshold from which a parameter gradient is considered significant. It can therefore be used to differentiate between a smooth parameter distribution and a sharp contrast that should be enhanced by the algorithm. Increasing $a$ decreases the influence of small parameter contrasts on the corresponding weight. For increasing roughness, the weighting becomes increasingly small reaching a minimum of $b$. Therefore, $b$ can be used to control the maximum increase of the gradients in one iteration. Depending on the choice of $a$ and $b$, the parameter roughness $r$ has the following effects on the weighing:

•	Small values for $r$ lead to a weighting of 1, i.e., no effect on the next iteration.
•	High gradients at boundary $i$ lead to small weights by locally lowering the global smoothness constraint and thus reduce the penalty at this boundary.
•	The behavior of the weighting function between the limits $b$ and 1 controls the convergence of the parameter gradients throughout the inversion.

Hellman et al. (2017, Figure 2) demonstrate the influence of the coupling parameters on the weighting function. However, for inversion, they choose standard values of 0.1 for $a$ and $b$ and do not try to optimize those values. Note that Hellman et al. (2017) introduce a third coupling parameter $c$. They leave this exponential term at 1 in their investigations, and, therefore, we decided to neglect this term for reasons of simplicity.

Finally, to exchange structural information between the parameter distributions, for each parameter a combined weight is calculated on basis of all other parameters. This ensures convergence of gradients at points, where other parameters also detect gradients, while omitting self-reinforcment. For a total of $M$ parameters to be coupled, the new regularization operator of each parameter $n\in M$ is determined as follows:

Figure 1 shows the flowchart for the general SCCI scheme, and Figure 2 demonstrates the effect of the different steps on a very simple two-layer model using the standard values of $a=b=0.1$. The start of the SCCI is always an Occam inversion of the data for MRS and VES separately. We use the result of the VES inversion to calculate the kernel function for the MRS inversion. All three models resemble a smooth representation of the sharp interface in the synthetic example (Figure 2a). We now start exchanging information between the parameters through structural coupling. The Occam inversion result is used as starting point of the coupled inversion, whereas the model-roughness distributions (equation 8) are used to reweight the inversion (Figure 2b).

For the presented two-layer case, the model roughness of the parameter distributions exhibit a local extreme close to the real model boundary at 20 m depth. Figure 2c shows the calculated weighting functions (blue curves) using equation 8. The combined weights (orange curves) are calculated according to equation 10 and used as local adjustment of the global smoothness constraint in equation 6 for the next iteration. As a result, the gradients start to increase once the coupling is started (Figure 2d, orange line). After some iterations, all three models show a sharp interface at the same position as a result of the coinciding model gradients due to the mutual weighting (Figure 2d, green line). Besides the sharp interface, the other parts are observed to be significantly smoother.

SYNTHETIC STUDIES

After the first simple case used to explain the principles of coupling, two more realistic models are used to outline the possibilities of structural coupling for different hydrologic environments. For both synthetic models, the influence of the parameters $a$ and $b$ are evaluated by the root-mean-squares (rms) between the SCCI result and the synthetic model for all parameter combinations $a$ and $b$. We vary each parameter in a logarithmic range between 0.01 and 1 with 10 steps per decade. For resistivity and relaxation times, a relative rms is calculated, whereas for water content, an absolute rms is used. Some extracted SCCI results are discussed in detail to give an idea of the effects of the coupling parameters with respect to the synthetic model behavior.

Case 1: Freshwater, clay, and saltwater

Whether used for exploration or monitoring, any geoelectrical sounding (VES, TEM, or HEM) has issues with distinguishing clay layers and saltwater aquifers. The first model consists of a unsaturated zone, a freshwater aquifer, a clay layer, and a saltwater aquifer beneath. It is a case typical for investigations in more shallow hydrogeologic environments at the coast and interesting for combining electrical and MRS because neither method is capable of explaining the whole model on its own. VES has difficulties to determine the difference between the saltwater aquifer and the clay, whereas MRS is unable to distinguish between saltwater and freshwater (Table 1). The used parameter set for the synthetic model is given in Table 2. The synthetic model is purely blocky with no smooth parameter contrasts.

For the measurement setup, we assumed 50 AB/2 spreads ranging logarithmically from 0.1 to 1000 m for the VES, and a 50 m diameter coincident circular loop with 20 pulse moments ranging from 0.1 to 10 As for the MRS. A relative amount of 3% Gaussian noise is added to the VES data, whereas 40 nV Gaussian noise is added to the ungated MRS data cube (the real and imaginary parts). The simulation and inversions are done on different meshes. For the inversion of the layered earth case, we used 47 layers to a depth of 50 m with a minimum thickness of 0.5 m for the first layer. The thicknesses increases following a hyperbolic sine function as proposed by Behroozmand et al. (2012b).

Figure 3a–3c shows the resulting rms values for all three inverted parameters. Because there is no MRS signal from the clay layer due to technical limitations, the layer is omitted in the rms calculation. The parameters $a=1$ and $b=1$ are used as a reference because this combination does not influence the coupling effect and therefore leads to the smooth model. We observed local minima of the rms (i.e., models closer to the true synthetic model) to be decreased by a factor of $2-4$, compared with the smooth inversion results. The rms is decreasing with decreasing $a$ or $b$ up to a certain point, where the parameter $a$ is too low. Lower values for $a$ can also lead to good rms values if $b$ is chosen high enough, especially for the MRS results. The lowest rms values are reached for $a$ between 0.1 and 0.5 for all three inverted parameters.

As the observed variation in the rms tends to be less dependent on the choice of $b$, we decide to take a closer look on the effect of parameter $a$. We use the summed quadratic roughnesses (total roughness) of the final SCCI results, over the coupling parameter $a$ for a low value of $b=0.025$ (Figure 3d), as additional tool for evaluating the different results. The summed values show how blocky our models are, with the value for $a=1$ being closest to the smooth model. We can observe a local maximum for values of $a$ approximately 0.3 and a local minimum for values close to 0.05.

We discuss the results of the SCCI more detailed for $a=0.215$, 0.046, and 0.022, marked with white crosses in Figure 3a–3c or black dots in Figure 3d. Figure 4 shows the SCCI results of the three different coupling parameter settings (from gray to black with the increasing number of iterations), the Occam type inversion (orange) and the synthetic model (blue). In comparison to the VES, the initial smooth inversion of the MRS is able to resolve the four layers of the synthetic model, but none of the models can be used to provide reliable information concerning the thickness and parameter values of the different layers. The SCCI result in Figure 4a with $a=0.215$ corresponds to the observed maximum in the roughness plot (Figure 3d). Therefore, it is the most blocky model found by SCCI. The result is very close to the synthetic model for all three parameters, with the exception of MRS, which is not able to resolve the clay layer in terms of parameter magnitude. This is an expected behavior because the water content of the clay layers is not detectable due to the very low relaxation times. Nevertheless, the layer boundaries are resolved quite well. We observe the lower boundary of the clay to be less well-resolved compared with the upper boundary. In contrast, the second SCCI result (Figure 4b), corresponding to the local minimum in the roughness plot, shows a very smooth behavior. The differences to the smooth inversion results occur at the synthetic-model boundaries, although we observe this model to be less affected by the structural coupling. Both results have an improved rms value compared with the smooth inversion result. The third model (Figure 4c) with $a=0.022$ has a higher rms than the original smooth inversion. It shows a high variability in all parameter distributions, especially for the water content. Compared with the other two results, the parameter estimation is difficult with a distribution similar to the smooth inverison, for constaints that are not chosen high enough.

Case 2: Vadose zone

In the previous case, we neglected a potentially smooth behavior of a vadose zone. Modeling a subsurface with smooth transitions (such as the vadose zone) and blocky features (layers) beneath is difficult with the block and Occam inversions. However, many field cases actually can only be described accurately if a smooth behavior is taken into account. Thus, our second case presents a more shallow investigation of the vadose zone. This synthetic example shows the effects of structural coupling if no major blocky features are present in the data.

For a realistic synthetic model, we calculate a water-retention curve after the scheme of Brooks and Corey (1964). They describe the effective saturation $(S)$ of a capillary transition zone over fully saturated aquifer as a function of the height $h$ by

with ${\lambda }_{BC}$ being the pore-distribution index and $h_0$ being the thickness of the fully saturated zone above the aquifer due to capillary forces. The resulting saturation based on the residual water content ${\theta }_R$ and the water content at saturation ${\theta }_S$ is then translated into a water-content distribution above the aquifer via We use a residual water content of ${\theta }_R=5\%$, a capillary elevation of $h_0=12\mathrm{cm}$, a saturated water content (porosity $\mathrm{\Phi }$) of ${\theta }_S=30\%$, and a ${\lambda }_{BC}$ value of 1; the water table is assumed at a depth of 2 m.

The electrical-resistivity distribution $\rho$ is calculated according to the second Archie equation (Schön, 2015) under consideration of the saturation $S$:

We chose a cementation exponent of $m=1.3$ (unconsolidated sand) and the saturation exponent of $n=2.0$. The water resistivity ${\rho }_w$ is set to $20\Omega \mathrm{m}$.

Costabel and Yaramanci (2011) provide a relation between relaxation time and saturation:

We chose a relaxation time of $T_2^{*}=200\mathrm{ms}$ for full saturation. Again, a relative amount of 3% Gaussian noise is added to the VES data, where in this case 2 nV Gaussian-distributed noise is added to the ungated MRS data cube (real and imaginary parts), leading to a signal-to-noise ratio of approximately 8. Again, the inversion and simulation are done on separate discretizations. For the inversion, we use 43 layers to a depth of 15 m with a minimum thickness of 0.1 m.

Similar to the layered case, Figure 5 shows the rms of the relative (resistivity, relaxation times) or absolute (water content) deviations as a function of $a$ and $b$. For $b$ between 0.01 and 0.1, the effects of the structural coupling are similar to that observed in the layered case. However, the optimum choice for $a$ (i.e., the minimum in the rms value) is shifted to smaller values, which increases the influence of small parameter gradients in the constraint weight. Because the vadose zone naturally consists of many small gradients, the structural coupling produces better results. In contrast to the layered case, another minimum for the rms values can be observed using high values of $b$ and smaller values of $a$. This leads to models close to that of the smooth inversion result, which already explains the smooth behavior of the synthetic model. In contrast to the layered case, the increase in rms values observed for low values of $a$ is smaller. Again, we decide to settle with a low value for $b$ and calculate the total roughness over parameter $a$ (Figure 5d). Similar to the previous case, a local maximum of the accumulated roughnesses can be observed, but this time for values of $a$ between 0.2 and 0.3. In contrast to the other case, the values do not concur with the minimum rms values in Figure 5a–5c.

Three representative SCCI results are discussed and presented in Figure 6 (from gray to black with the increasing number of iterations). The smooth inversion result (orange) and the synthetic model (blue) are shown for comparison. Figure 6a represents the results of the SCCI for $a=0.063$ and $b=0.025$. The parameter trend resembles the synthetic model and leads to two to three times smaller rms values compared with the smooth result. The Occam-type inversion is able to find the correct parameter beneath the vadose zone, but it has difficulties resolving the correct parameter trend in the first 2 m. The curvature of the Brooks-Corey transition zone is fitted much better by the structurally coupled inversion. Also, the SCCI result reaches the parameter values for the saturated aquifer earlier than the smooth result. Additionally, the sharp boundary at the water table is resolved more clearly. The second result corresponds to $a=0.292$ (Figure 6b) and is observed to be the local maximum in the presented roughness plot (Figure 5d). It shows the highest parameter contrast of all SCCI results and can be compared with the results of a Marquardt inversion. Still, this result has a slightly better rms compared with the smooth inversion. Figure 6c shows the results within the second mentioned minimum in the rms plots, corresponding to a high value of $b=0.341$ and $a=0.063$. Compared with the other results, the model update from one iteration to the next is small. Therefore, the SCCI result is close to the original smooth inversion result. The parameter estimation of the resistivity and relaxation times are similar to that of the model presented in Figure 6a, except for the sharp bend, at which the saturated aquifer begins. However, also similar to the smooth result, the water content is overestimated for the transition zone.

Investigating the smooth vadose zone in the first 3 m yields that a small value of $a<0.1$ leads to better results than higher values. Note that all models, regardless of their rms, are valid as long as they have an equivalent data fit compared with the smooth inversion.

FIELD EXAMPLE

A well-understood field data set from the island of Borkum, Germany, has been chosen for field-scale validation. The data were previously published in Günther and Müller-Petke (2012) and were inverted using a joint-inversion approach with a blocky discretization. From all presented soundings in the original paper, we choose the first (CL2) because it shows the most complicated subsurface and the worst (most realistic) data quality. Ground truth is available from the research drillhole CLIWAT-2 and array induction logs (Figure 7). It clearly shows a smooth transition from freshwater to saltwater greater than a 20 m depth interval for which a block discretization is inappropriate, although this is possibly beyond resolution limits. It becomes clear that the presented five-layer model cannot adequately describe the geology in sufficient detail. Particularly, the middle silt layer is not resolved as shown by the error bars in Günther and Müller-Petke (2012, Figure 3).

The MRS was done with a rectangular loop of $49\times 47\mathrm{m}$ size, using 36 pulse moments that are logarithmically distributed between 0.07 and 7 As. The noise level of the data has been reduced using a reference coil and noise-canceling techniques. Note that Günther and Müller-Petke (2012) compute a circular loop of equivalent area due to the nonavailability of a code for rectangular loops. Additionally, they do not use resistivity simultaneously for kernel calculation but only once at the beginning of the joint inversion and again toward the end. Because the data do not provide correct phase information, we inverted rotated amplitudes. The VES consists of 23 measurements with current half-spread AB/2 ranging logarithmically from 1.5 to 150 m. For the potential, half-spreads of MN/2 = 0.5 and 5 m are used.

The hydrogeologic setting is represented by a vadose zone and three aquifers, separated by several embedded aquitards. The first aquifer (well-sorted aeolic sands) has a thickness of 23 m with a water table at a depth of 3 m below the surface. Underneath is a silt layer of a few meters thickness followed by alternating sequences of fine sand and clay to a depth of 32 m. The second aquifer from 32 to 49 m also consists of (poorly sorted fluviatile) fine sands and contains fine clay layers in a band of a few meter thickness.

The single inversion of the VES indicates the beginning of the freshwater aquifer at a depth of 3 m. The profile shows a smooth transition zone, which cannot resolve the two clay layers, indicating resistivities of approximately $100\Omega \mathrm{m}$ at the beginning of the aquifer and values of $10\Omega \mathrm{m}$ at the bottom. The water content varies between 25% and 30% with a maximum at a depth of 15 m. In contrast to that, the expected water content from the drilling in the second aquifer can be assumed to be 30%–35% and the water content of the clay sand layer should drop dramatically because the water content of the clay is not detectable due to the low relaxation time. The relaxation time, however, shows very little contrasts varying between 200 and 300 ms.

Based on the observations of the synthetic models, a low value for $b=0.05$ is used while investigating the different results for a range of logarithmically distributed values for the coupling parameter $a$. Lacking a synthetic model for comparison, we rely on the total quadratic roughness of the combined models to evaluate the SCCI results (Figure 8). Similar to the synthetic modelings, high values of $a$ have no influence on the inversion result. Starting coupling in the SCCI leads to strong contrasts in the water content and relaxation times in the first meters as well as at the basis of the first aquifer. Reducing the values of $a$ quickly leads to blocky models. The model corresponding to the first local maximum of the total roughness $(a=0.089)$ is plotted in Figure 7 (light blue). The model is very close to the block inversion result of Günther and Müller-Petke (2012).

Further lowering of $a$ increases the number of gradients supported by the structural coupling leading to smoother results, especially for the water content. As a second result, we show the model for $a=0.01$ (Figure 7, dark blue line). Similar to the results of the synthetic modeling, a low value for $a$ leads to more variation of the model. Without further knowledge, this model would possibly be discarded in favor of the other (smoother) model. However, the presented features actually coincide much better with the lithology than the blocky result. We can observe lower water contents and relaxation times in the aquitard region (23−32 m), but also an increase of water content at the bottom of the upper aquifer (coarser material) and a decrease of water content toward the bottom of the second aquifer (clay layers reducing integrated water content).

Note that all models show a possible solution for the presented field case. Changing the value for $b$ to 0.025 or 0.2 results in similar features, even if the presented figures and discussion are restricted to models with $b=0.05$. We decided to show extreme models: the most blocky and the most smooth. The real model has to be included in this range, as shown with the synthetic models. Both models agree with the additional information and are able to resolve the second aquifer and the first silt layer in terms of water content and relaxation times. The resistivity sounding cannot resolve the first clay layer in 25 m depth if the model is too blocky. The base of the second aquifer at a depth of 50 m is detected in all three parameter distributions and in all models between the blocky and the smooth solution. Keeping in mind the loop layout, a depth of 50 m is not as well-resolved as the first contrast in 20 m. Compared with the initial smooth inversion results, the coupled models align better with the expectations from the borehole lithology. The water content is concentrated in the aquifer, rising to values up to 30%, whereas at less than 50 m, the values are reduced compared with the smooth result. The major contrasts due to changes in the lithology are resolved a little more clearly in the blocky model. To sum it up, particularly due to the improved water-content estimation in the deeper part, all layers identified in the drilling can be found in the SCCI result.

Figures 9 and 10 show the misfit of the two presented SCCI results for MRS and VES, respectively. The remaining misfit of the blocky and the smoothly varying model only contains Gaussian distributed noise. The observed misfit is representative for all models shown and not shown in this study, including those from the synthetic studies.

DISCUSSION

The results of the SCCI range between the models that a block or smooth inversion can produce, while providing a similar data fit. However, the best fit in terms of rms compared with the synthetic model is neither the most smooth nor the most blocky solution. The SCCI therefore is able to find models that are much more consistent with the synthetic models, or lithology. For a subsequent block inversion, the results can be used to determine the number of needed layers and to generate starting models. In the case of the vadose zone, the SCCI can be used to fit water-retention parameters directly.

The exchange of structural information does not only affect the parameter distributions directly. Implicitly, the improved resistivity distribution changes the kernel function due to the magnetic-field calculation. It stands to reason that this is beneficial for complex MRS inversion, whose dependency on an accurate resistivity distribution is greater compared with rotated amplitudes (Braun and Yaramanci, 2008).

Effects of the weighting parameters

The coupling is controlled by two coupling parameters. We show that there are main parameter contrasts that can be observed throughout all the different models and some features are only seen in models for specific combinations of $a$ and $b$. Different coupling parameters correspond to either mostly smooth or mostly blocky models. Similar to a minimum-gradient-support inversion (Grombacher et al., 2017), a smooth parameter distribution can be resolved together with sharp contrasts, which is neither possible for block nor for smooth inversions. The exchange of structural information can lead to more reliable estimations of layer depth and thickness even if a single method cannot resolve a layer boundary in their correct shape. This can be observed in the synthetic model of the vadose zone, in which both methods tend to see the beginning of the saturated aquifer in greater depth than the synthetic model.

Generally, higher values for $a$ lead to a constraint weight that is only affected by the major gradients in the smooth starting model. Thus, for high values of $a$, higher parameter contrasts are needed to affect the coupling weight. Values for $a>1$ do not lead to any kind of coupling, although the value depends on the roughness distribution which is specific for each geophysical method (and model transformation). Reducing $a$ at first leads to more blocky models resulting in a local maximum of the quadratic roughnesses, seen in the layered and the vadose-zone case. Throughout the SCCI iterations, the major gradients are able to focus on one boundary instead of being smeared over many. Note that at this point, we did not investigate the number of iterations but we expect that similar results can be obtained with different $a$ and $b$. Concerning the example of the vadose zone, the contrast between full and partly saturated material is the only contrast affecting the coupling, leading to models with only one main boundary, which is the water table. A further decrease of $a$ leads to small roughnesses because more parameter contrasts are affected by the coupling. This usually leads to models close to the smooth inversion result, corresponding to a local minimum of the quadratic roughnesses. In between those two extrema, all models between blocky and smooth can be found. Depending on the roughness of the initial smooth inversion result, very small values for $a$ do not lead to reasonable results even if the data fit is satisfied. For example, for $a=0.022$ (the layered case), nearly all gradients in the model have significant influence on the combined constraint weight. The overall reduction of the smoothness constraint destabilizes the inversion and results in highly oscillating models. Note that all SCCI results are (such as the smooth inversion) valid results of the inverse problem, and models with a better rms are actually improved compared with smooth inversions.

Parameter $b$ stabilizes the inversion in a certain range of parameter $a$. For $1\ge b\ge 0.5$, only strong parameter contrasts have an influence on the constraint weight and the model update for each iteration is limited. The initial smooth result remains nearly unchanged. For $0.1\le b<0.5$, a general shift in the optimum value of parameter $a$ to smaller values can be observed. Best results are obtained for $b<0.1$ to stabilize the inversion and searching for the optimal $a$.

In contrast to the layered case, the vadose-zone case is observed to have a second maximum for high values of $b$, which does not behave like a simple shift of the correct $a$ as mentioned above. The minimum corresponds to low values of $a$ and high values of $b$. Therefore, the results are close to the smooth initial result, which is, in the case of the vadose zone, already a very good model for the synthetic case. Supporting most gradients leads to slightly stronger gradients at the bottom of the transition zone, hence the better rms. However, the range of appropriate coupling parameters is narrow and might be difficult to find if no synthetic model (or other methods for validation) is at hand.

Because structural coupling locally reduces the global smoothness constraint, a slight reduction in ${\chi }^2$ values can be observed. This is mainly caused by the first iteration step of the structural coupling, in which the strongest reduction of the constraint due to the weighting can be expected. This is also the reason why very small values of $a$ lead to inappropriate models. The oscillating models are overfitting the data, resulting in ${\chi }^2\ll 1$. In this case, higher values for $a$ lead to better results. After the first coupled iteration, the model roughness becomes more focused on single boundaries over the iterations and the overall reduction of the global smoothness constraint becomes negligible. Nevertheless, the data fits of all models only contain white noise and therefore are able to fit the data within the overall noise level.

The use of a combined weight can lead to very small values in the coupling for later iterations when coinciding boundaries are found. A further decrease of the regularization is not needed at that point of the inversion. To avoid instability, we set a lower boundary of 0.04 for the combined weights. The applied boundaries for the minimum and maximum combined constraint weights (0.04 and 1, respectively) are chosen independently of the coupling parameters $a$ and $b$ and lead to a maximum factor of 25 between the lowest and highest constraint at each iteration. The maximum value results in the same penalty for a parameter contrast as in the smooth inversion. During the parameter study, the minimum value only comes into account in inversions using the lowest coupling parameters or late iterations and therefore does not alter the results nor conclusions derived by the study.

Recommendations for choosing coupling parameters

As a first interpretation for choosing the best parameter combination for a case, we observe that the range for a reasonable $a$ is smaller than the range for appropriate $b$. However, while choosing $b>0.1$, not all features of the synthetic models are resolved because the overall change from the smooth model per iteration is limited. Generally, we recommend setting $b$ to a small value (<0.05) to stabilize the inversion and then search for a value for $a$ to find a geologically reasonable model. The summed quadratic roughness can be used to identify the most blocky and the smoothest inversion result, to get an idea of the model. That can be seen at the layered-case model, in which the summed quadratic roughness identifies the most blocky model, with the lowest overall weighting, similar to choosing the correct parametrization of a block inversion. The variability of the models as a result of the SCCI contains all models between block and smooth, as shown in the field case. The models contain all information from the single data sets as well as the additional structural information. All models produced with the structural coupling are valid in terms of data fit, leaving it to additional information to find the best model or coupling parameter for the actual problem at hand.

The optimum value for the coupling parameter $a$ is shifted to higher values, if a lot of sharp boundaries are observed in the model. For $b\approx 0.05$, a broad range of $a$ values leads to improved models. With increasing values of $b$, the range for the best $a$ becomes narrower and therefore more difficult to find.

Insight into the space of equivalent models between smooth and blocky characteristics can be achieved by using a range of SCCI parameters. We recommend a fixed value for parameter $b$ and a logarithmic distribution for parameter $a$. The total roughness can be used to choose some representative models for the final hydrogeophysical interpretation.

CONCLUSION

We applied the concept of SCCI to the joint inversion of electrical and MRS using a fixed model discretization with smoothness constraints. Compared with individual inversions, the exchange of structural information between resistivity, water content, and relaxation time leads to more contrasts in the models if the parameters show coinciding gradients. The algorithm can be applied to blocky or smooth parameter variations or models showing smooth transitions and sharp changes. The coupling is controlled by two parameters that can be used to incorporate the expectations on the parameter distribution. Whereas one can be fixed, the other should be varied to produce models of different complexity. We showed that inversion results of two synthetic models, a layered and a smooth case, are improved by using structural coupling. Moreover, most $a$-, $b$-combinations improve the parameter boundary estimation of the smooth inversion and lead to models closer to the synthetic model. Additionally, layer boundaries are more distinct while the parameter distributions, apart from big contrasts, are smoother than in the starting smooth inversion. In the case of a layered earth model, all parameter distributions show consistent boundaries for the top and the bottom of the first aquifer as well as the bottom of the clay. Structural coupling does not force parameter contrasts that are not supported by the data. However, they can occur if the parameters are not chosen carefully, for example, for very low values for $a$. We have shown that SCCI can effectively be used to reduce the ambiguity of the different methods. As a result, one obtains all equivalent models between smooth and blocky that fit the data, which gives a better and more objective understanding of the investigated target. Even if we only invert 1D VES and MRS problems, the presented method for coupling can easily be applied to other types of soundings, but also extended to a 2D structurally coupled inversion.

DATA AND MATERIALS AVAILABILITY

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