## ABSTRACT

Nuclear-magnetic resonance (NMR) is a powerful tool for groundwater system imaging. Ongoing developments in surface NMR, for example, multichannel devices, allow for investigations of increasingly complex subsurface structures. However, with the growing complexity of field cases, the availability of appropriate software to accomplish the in-depth data analysis becomes limited. The open-source Python toolbox coupled magnetic resonance and electrical resistivity tomography (COMET) provides the community with a software for modeling and inversion of complex surface NMR data. COMET allows the NMR parameters’ water content and relaxation time to vary in one dimension or two dimensions and accounts for arbitrary electrical resistivity distributions. It offers a wide range of classes and functions to use the software via scripts without in-depth programming knowledge. We validated COMET to existing software for a simple 1D example. We discovered the potential of COMET by a complex 2D case, showing 2D inversions using different approximations for the resistivity, including a smooth distribution from electrical resistivity tomography (ERT). The use of ERT-based resistivity results in similar water content and relaxation time images compared with using the original synthetic block resistivity. We find that complex inversion may indicate incorrect resistivity by non-Gaussian data misfits, whereas amplitude inversion shows well-fitted data, but leading to erroneous NMR models.

## INTRODUCTION

The surface nuclear-magnetic resonance (surface NMR) method has proven to be a valuable technique in hydrogeophysical investigations thanks to its direct sensitivity to the hydrogen proton and, thus, subsurface water. Most surface NMR applications target the exploration of 1D structures; that is, they deliver depth distributions of the derived NMR parameter water content and relaxation time. Various case studies can be found in the literature, including permafrost investigations (e.g., Parsekian et al., 2013; Keating et al., 2018), vadose zone characterization (e.g., Roy and Lubczynski, 2005; Costabel and Günther, 2014), karst environments (e.g., Vouillamoz et al., 2003; Mazzilli et al., 2016), or salt-water intrusion (e.g., Günther and Müller-Petke, 2012; Costabel et al., 2017). In recent years, however, surface NMR has been developed toward the investigation of 2D and more complex structures. In this context, the technique is commonly referred to as magnetic resonance tomography (MRT). The development of MRT started with a sensitivity and resolution study introducing the separated transmitter and receiver loop approach by Hertrich et al. (2005), based on the general surface NMR forward formulation presented by Weichman et al. (2000). Next, inversion schemes and field data examples were presented (Hertrich et al., 2007, 2009; Dlugosch et al., 2014). The development of MRT has been supported by the availability of multichannel measurement devices (Walsh, 2008; Dlugosch et al., 2011; Liu et al., 2019). Dedicated 2D survey layouts such as the elongated transmitter and in-loop receiver array (ETRA) after Jiang et al. (2015b) allow for reasonable measurement times to be achieved. However, even though 2D NMR structures are imaged, the subsurface resistivity is mostly modeled as a 1D structure. Hence, another step toward a comprehensive 2D investigation is the possibility to calculate magnetic fields on the basis of 2D resistivity distributions. The influence of 2D resistivity subsurface structures on surface NMR measurement has been the subject of discussion in several publications (e.g., Braun and Yaramanci, 2011; Lehmann-Horn et al., 2011, 2012). Common to all studies are block-like resistivity structures, where several subsurface blocks are defined and a single resistivity value is assigned to each block. Whereas using block-like resistivity structures appears reasonable for synthetic studies, in field cases, mostly smooth resistivity distribution is obtained by electrical resistivity methods; thus, an additional step generating block-like elements is necessary. Consequently, it is favorable to directly include smooth 2D or even 3D resistivity structures directly into the MRT forward calculation, to avoid this step.

Besides 2D MRT inversion, the joint inversion of layered 1D surface NMR and resistivity provide superior resolution and information (Günther and Müller-Petke, 2012; Behroozmand et al., 2013). Skibbe et al. (2018) show that a coupled surface NMR and resistivity inversion based on a smooth distribution is more flexible as it allows for the inclusion of structures such as the capillary fringe or salt-/freshwater mixing zones. In this context and when targeting 2D joint or coupled inversion of MRT and resistivity, again, it is favorable to be able to handle smooth resistivity models.

The topic of reproducible science has gained attention in recent years (e.g., Broggini et al., 2017). The idea is to ease the way for others to retrace all the steps from the original data to final results. Besides revealing the underlying data, the availability of the used software is essential for sustainable research. Open-source software has made its way in geophysics, particularly, in electromagnetics. For geophysical applications, interpreted languages such as MATLAB (e.g., Key, 2012), Python (e.g., Werthmüller, 2017), and Julia (e.g., Han et al., 2018) are very common. Python’s popularity in the geophysical community has resulted in a large ecosystem of open-source packages and tools (e.g., Cockett et al., 2015; Rücker et al., 2017).

Concerning surface NMR there is a range of available and in particular open-source software. To the best of our knowledge, there are three available codes.

- •
The AarhusInv package (Auken et al., 2015) includes a surface NMR package (Behroozmand et al., 2012), providing data processing, modeling, and inversion in one dimension. The software is free for academic use but not open-source.

- •
The Akvo (Irons et al., 2019) package interfacing with Merlin a surface NMR module that provides data processing, 1D modeling, and 1D inversion. Sources are available at GitHub.

- •
MRSmatlab (Müller-Petke et al., 2016) provides data processing, modeling, and inversion in one dimension. Sources are available at GitHub.

Furthermore, there are proprietary software packages provided by the instrument manufacturers (Vista Clara and Iris Instruments) that are delivered with the instrument but are not open source. All of these codes handle 1D problems; some of them handle 2D problems as well. However, algorithms for 2D investigations are presented in several publications, but open modeling and inversion code are not freely available. Beyond that, there is no publication that accounts for more complex resistivity structures than fixed 2D blocks.

In this paper, we introduce coupled magnetic resonance and electrical resistivity tomography (COMET), the first open-source toolbox capable of MRT forward and inverse problems in two dimensions with arbitrary resistivity structures. COMET also includes the routines already presented in Skibbe et al. (2018) and thus provides 1D capabilities as well.

## SOFTWARE DESIGN

COMET is developed in Python and supports versions 3.6 and newer, as the used third-party packages. Python allows us to gather functionalities and input data in different scopes (modules and classes). From these central places, information is shared for the surface NMR modeling and magnetic field calculation. This avoids redundancy and simplifies the more complex tasks. COMET includes 1D and 2D capabilities for surface NMR modeling and inversion and is designed to be easily extended to three dimensions. COMET does not include data processing (aside from gating) because there are tools well-suited for this purpose, such as MRSmatlab (Müller-Petke et al., 2016). Import routines for MRSmatlab data files are provided, but COMET can deal with any kind of numerical arrays, for example, those provided as simple ASCII files or the like. Running COMET requires, apart from common libraries included in major python distributions, two third-party packages: pyGIMLi (Rücker et al., 2017) and custEM (Rochlitz et al., 2019).

The pyGIMLi package is the basic framework that we used to develop COMET. We use a broad spectrum of methods from pyGIMLi including but not limited to mesh and marker handling, matrix algebra, as well as the Gauss-Newton algorithm with incomplete line search that we use within the inversion routine. With the latter, COMET is able to invert water content and relaxation time distributions using a Levenberg-Marquardt scheme (block inversion) in one dimension or a smoothness-constrained inversion in one dimension or two dimensions. Skibbe et al. (2018) introduce a structurally coupled inversion approach between water content, relaxation times, and resistivity in one dimension using the COMET toolbox. Whereas the general-purpose library pyGIMLi is always required, custEM is needed only in cases in which the 3D magnetic field computations build upon a non-1D resistivity distribution. The custEM package is built upon the open-source finite-element library FEniCS (Alnæs et al., 2015). As a default mesh builder in two dimensions, we use Triangle (Shewchuk, 1996). In three dimensions we rely on TetGen (Si, 2015) as the default. However, external meshes can easily be imported and used without the requirement of a TetGen installation.

The main features of COMET are gathered in class objects. These are the main interfaces for the user, connecting the input parameters with the automatic routines that handle most of the processing. Figure 1 introduces the main features of COMET and possible user interactions from the source discretization to the inversion results.

There are four classes as main interfaces users deal with. Additional calls from other parts of the COMET source code can usually be avoided, which makes COMET easy to use, even if the user has little Python experience. We divided COMET in two independent parts (apart from their organization in software classes): first the magnetic field calculation of the transmitter (tx) and receiver (rx) loops and, second, the forward calculation of the NMR signal. These two parts are described in the following.

### Magnetic fields

The first submodule handles source configurations and the magnetic field calculation. COMET provides instances of loop classes to handle different source types. The loop class represents the source configuration and contains all functionalities for the magnetic field calculation including but not limited to receiver/transmitter discretization import and creation as well as resistivity handling in one dimension or two dimensions (Figure 1). In detail, the module discretizes any surface loop with horizontal electric dipoles (HEDs) and calculates the summed magnetic (or electric) field of the dipole combination to describe arbitrarily shaped loops segment-wise. The basic formulations used for the magnetic field calculations have been published numerous times, going back to Ward and Hohmann (1988) or more recent implementations of Key (2009) or Werthmüller (2017). The calculation of magnetic fields for 1D-layered media is already well-investigated and is therefore not the focus of this manuscript.

A pure Python implementation for primary electromagnetic fields gives COMET several advantages because it allows complete 1D modeling and inversion with just one package. Furthermore, it provides an easy and robust meshing of arbitrary source types with respect to the requirements of a surface NMR calculation, for example, in the vicinity of the loop wires. Each used loop (tx and/or rx) is represented by an instance of the loop class, which comes with numerous methods for adapting the sources and source discretizations. This means that COMET can calculate the fields on any discretization the user chooses or it can automatically generate a suitable discretization for the purpose of surface NMR. These instances are then distributed and used by other parts of COMET. Also, the mesh and corresponding fields are stored by the loop classes and can be saved and loaded for later usage.

However, if the resistivity distribution is more complex than one dimension, COMET provides those magnetic fields as primary field and uses the secondary field approach from the third-party module custEM to calculate the total field. Rochlitz et al. (2019) introduce several approaches to solve the Maxwell equations, but as a default, we rely on the magnetic field approach with the advantage of omitting the calculation of electric fields. This saves time with respect to other approaches typically applied in EM modeling. Furthermore, this part of COMET is used as the primary magnetic field calculation routine in custEM and validated against other codes by Rochlitz et al. (2019).

### Surface NMR

We use the surface NMR basic equations (e.g., Weichman et al., 2000; Hertrich, 2008) to calculate the kernel function. For separated loops, we implement the formulations of Hertrich et al. (2005). To account for off-resonance effects, we use the formulations of Walbrecker et al. (2011). COMET is able to invert single coincident or separated loop soundings or any combination of soundings jointly, but it is currently limited to single-pulse free-induction decay (FID) measurements. Water content and relaxation times are inverted according to the qt-scheme introduced by Mueller-Petke and Yaramanci (2010) that inverts the whole data cube spanned by pulse moment and measuring time simultaneously. However, we constrain the relaxation time distribution to be a monoexponential decay (smooth-mono or block-mono) in every model cell (Dlugosch et al., 2014; Müller-Petke et al., 2016). COMET supports the complex inversion of surface NMR data, that is, real and imaginary parts of the data are fitted independently, in the time domain (Weichman et al., 2000; Braun and Yaramanci, 2008). However, this approach can lead to difficulties if the phase of the NMR data cannot be fitted (Grombacher et al., 2016; Grombacher and Auken, 2018). In those cases, COMET offers the possibility to invert rotated amplitudes instead. Therefore, the data set is rotated until the power of the imaginary part is minimized and the real part of the data is inverted via an amplitude inversion. Mueller-Petke and Yaramanci (2010, Appendix B) state the necessary alterations of the Jacobian matrix for amplitude inversion. Müller-Petke et al. (2016, Appendix D) and Jiang et al. (2017) give the Jacobian matrix for complex inversion in one dimension and two dimensions, respectively.

Due to the complexity of tasks, three different classes are used to handle the surface NMR modeling (Figure 1). The main object to deal with is the survey class, as it contains all information about sources, measurement configurations, and the data sets itself. The survey class can either be created bare and filled with all needed information (e.g., using native lists and arrays) or is created automatically importing an appropriate file (e.g., an mrsd file from MRSmatlab or a previously saved survey). Most other classes only need the survey class as the basic input for further calculations and use it to store everything from a single measurement to a whole complex 2D setup. Currently, COMET supports only the FID that is represented by an FID class instance. The survey class holds references to the (multiple) loop and FID classes and connects them.

The second important task is handled by the kernel class, which is mostly created automatically by the FID classes. A kernel class needs a reference to a survey and a simple index to the FID it represents. Through the FID, the kernel gets the needed information about the transmitter (and optionally a separate receiver) plus the corresponding magnetic fields as well as other information needed to calculate the kernel matrix, that is, pulse moments, frequency offsets, dead times, and others. The kernel class also provides methods for the automatic generation of suitable kernel discretizations in one dimension or two dimensions.

Manager classes for modeling and inversion are written in the same style as those of pyGIMLi to use them with the pyGIMLi inversion engine. We provide several derivatives of the basic manager class for specific purposes, for example, one dimension, two dimensions, or joint inversion (generalized as modeling in Figure 1). The main tasks are the forward modeling with a given kernel matrix and model and the computation of the Jacobian matrix for the inversion. The survey class provides the manager with the needed references to a corresponding kernel class per FID. Like the kernel class, the manager instances are created automatically as all the input parameters are already gathered within the loop, survey, and FID classes.

We leave the in-depth descriptions of the class interactions as well as specifics of the syntax to examples and tutorials provided with the source code. However, an example code that covers a simple 1D study is given in Appendix A.

## 1D VALIDATION

Because most parts of the COMET source code are written from scratch, they need to be validated using known examples and external software. Behroozmand et al. (2012, Figure 5) show a comparison of available codes including MRSmatlab based on a simple test case. The authors show that their code (AarhusInv) is in good agreement (deviations below 5%, mostly within 2%) with MRSmatlab. To evaluate COMET, we therefore select the same model, calculate a sounding curve, and compare it with the result of MRSmatlab (Figure 2). The synthetic model consists of three layers, the first two layers being 10 and 15 m thick with resistivities of 50 and $200\u2009\u2009\Omega m$ over a half-space with $20\u2009\u2009\Omega m$. The temperature is fixed to 293 K. A $100\xd7100\u2009\u2009m$ coincident loop is used to calculate the sounding curve for a 100% water content half-space. MRSmatlab uses an analytic expression for a circular loop with an equivalent surface area of $10,000\u2009\u2009m2$ and cannot calculate square loops directly. Both codes use the same $z$ discretization, but COMET has a different horizontal discretization because it uses 0.5 m dipoles for the loop discretization. We show the two sounding curves calculated with COMET for the square loop as well as for the circular loop and compare them to MRSmatlab.

Figure 2 shows that both codes produce very similar sounding curves. For the higher pulse moments, the absolute signal reaches a maximum of approximately 3500–4500 nV. This translates into a maximum relative deviation of approximately 2%–3% (100–150 nV). Taking into account that the sounding curves are calculated for 100% water content, the difference for a realistic case is less than the low-noise condition ($0.1\u2009\u2009nV/m2$) and well within the range of errors shown by Behroozmand et al. (2012). We assign the observed differences to the different discretizations.

## 2D EXAMPLE

One of the key targets of the COMET development is the 2D modeling and inversion of MRT including arbitrary 2D resistivity distributions such as output data from ERT inversion software. Note that COMET handles the parameter dimension between surface NMR modeling and inversion independent from the resistivity dimension; that is, a 2D MRT modeling can be based on 1D resistivity distributions and vice versa. This is because there are reasonable cases, such as salt-water intrusion into a coastal aquifer, in which either the surface NMR parameters are 1D, while the resistivity is 2D or where even a resistive half-space is sufficient, for example, for small loops. In the following, we present and illustrate the COMET capabilities for a complex 2D MRT. A special focus is set on demonstrating the impact of different resistivity implementations.

### Forward modeling

Figure 3 shows the synthetic model, whereas Table 1 summarizes the used parameters. We define two different aquifers interrupted by an aquitard. However, the first aquifer interrupts the aquitard at some point and connects to the second aquifer forming the structure of a buried valley. Furthermore, we define an unsaturated zone in the first 3 m. We explicitly avoided a symmetry with respect to the $x$ direction, so that the slopes of the buried valley show a different steepness.

With respect to the resistivity, there are different approaches to investigate such a target. We present and compare three of these approaches using the COMET toolbox. In all three cases, we invert for a 2D water content and relaxation time distribution, but we implement the resistivity (1) as blocky elements, (2) by a 1D simplification, or (3) using a smooth 2D distribution.

For the first case (the blocky case), we use the original synthetic parametrization. For this case, the resistivity is the same as for the synthetic data set; therefore, we can use this case as a reference to evaluate the other cases. The second case, a rough simplification constraining the resistivity to be one dimensional (the same values as before), ignores the 2D valley structure (the layered case). Third, we take a smooth distribution of the resistivity. This represents how to approach this problem under field conditions. For that, we simulate a synthetic ERT measurement with 121 electrodes in a 2 m spacing between $\u2212120$ and 120 m, using a multiple gradient measuring protocol. We use the implemented inversion routine based on pyGIMLi to invert our third input resistivity distribution for the MRT inversion shown in Figure 4 (the ERT case). We refrain from an in-depth discussion of the inversion or choice of parameters because ERT inversion techniques are well-investigated and are not a focus of this paper. The whole modeling workflow is illustrated in Figure 5.

Besides the resistivities, the MRT loop setup needs to be defined. We use an ETRA survey after Jiang et al. (2015b) with eight square loops of a base length of 50 m (Figure 5a). This means that our transmitter loop (black) is a $150\xd750\u2009\u2009m$ loop and we use seven additional 50 m square receivers (red and blue), distributed half overlapping in the tx loop. The edge length of 50 m provides the necessary resolution and investigation depth for the target structure.

Next, for calculating the kernel, we need the magnetic fields of the transmitter and receiver loops. In the example case, these are calculated with the finite-element method (FEM) on an appropriate mesh that is automatically created based on the input survey layout. Figure 5c shows part of the actual 3D tetrahedral FEM mesh with the transferred resistivities of the ERT inversion result in Figure 5b. In this case, we used 7812 different resistivity values in two dimensions, defined on triangles. For an accurate calculation, a typical FEM mesh needs to be larger than the 2D mesh. To avoid artifacts at the edge of the 2D mesh, we need to extend the resistivity values in the FEM mesh instead of letting the resistivity drop to the background model. The background resistivity is usually set to match the values near the source to have an accurate primary field in the vicinity of the loops. The values from the 2D mesh are transferred to the 3D cells using a nearest-neighbor interpolation. For one of the loops, the resulting total magnetic field obtained with custEM is shown in Figure 5d.

As a next step, we require the kernel function on an appropriate inversion mesh. Keeping in mind that the kernel is a highly alternating function in the vicinity of the source, usually many cells, typically more than 5 million nodes, are needed to obtain a reasonable representation. Together with the necessity for a clean integration in one dimension or two dimensions, a standard FEM mesh is rather unsuited for the direct calculation of a kernel function. However, the magnetic field, being continuous at resistivity changes, is easy to interpolate. Therefore, COMET does not calculate the kernel function on the same unstructured FEM mesh. Instead, we use a triangular prism mesh to ensure an appropriate integration. To accomplish this, an irregular 2D triangle mesh (in this case, the ERT inversion mesh) is built and then extended to three dimensions. This gives us complete control over the resistivity input and kernel output discretization without restricting the optimization of the FEM mesh (i.e., Nedelec elements are used in custEM), therefore the performance of our FEM solver. To obtain a reasonable representation of the kernel on a standard ERT inversion mesh, which is too coarse, we applied an internal refinement of the mesh cells (e.g., Jiang et al., 2015a). Here, we do not use numerical integration, but we summarize the kernel over refined subcells. More specifically, the triangular prism cells are refined by order $N$ so that each cell is divided in $4N$ prisms and the kernel function is estimated for every midpoint. For the final result, the refined kernels are integrated on the original cell discretization with an appropriate weighting. In the presented cases, we use an order of $N=2$. The kernel function itself can be calculated for each cell (or subcell) independently, which can be parallelized by distributing the 2D slices to different cores. Parallelization at this point reduces computation times as well as the memory consumption of the kernel calculation to a degree that it can be handled by any desktop computer.

The full 3D kernel, calculated for the cell centers of the 3D mesh using the interpolated magnetic field values from the FEM mesh, is shown in Figure 5e. Figure 5f illustrates the integrated 2D kernel for one tx-/rx combination and pulse moment. Note that the 3D kernel shown in Figure 5e is only visualized to understand the processing step. The kernel is instead calculated for the single 2D slices (for discrete $y$ values) in parallel and integrated (f) on the fly. Code examples for the ERT case are given in Appendix B.

This represents the most general case of the kernel calculation. For simpler cases, several abbreviations can be used. Calculating a 2D kernel for 1D resistivities is straightforward because the magnetic field calculation is simplified to a simple primary field calculation using semianalytic Hankel formulations. Calculating the kernel in one dimension is done on a similar triangle prism mesh, however with an extension in $z$ (using the inversion discretization). Hence, the integration is performed over the unstructured triangles in the $x$-$y$ plane such that the kernel is a function of $z$ only.

With the values in Table 1 and the kernel functions, the synthetic data are calculated using the true blocky resistivity model. We added 50 nV Gaussian noise to the separated loop data and three times as much to the coincident data of the transmitter loop due to the higher loop area. This data set is then used for all subsequent inversions.

### Inversion

To investigate the impact of the different resistivity approaches, three inversion kernels are calculated based on the three resistivity approximations, respectively. These are used to invert the synthetic data for water content and relaxation times using a rotated amplitude inversion as well as a complex inversion. Thus, we have six different inversion results and data fits (three approximations times two inversion strategies). For all cases, we use a Gauss-Newton inversion on the ERT inversion discretization. The regularization is controlled by a smoothness constraint $\lambda $, which is chosen in such a way that the remaining error-weighted misfit is Gaussian distributed with a root-mean-square of $\chi =1$. Additionally, an anisotropic regularization ($z$-weight) is used to control the smoothness in different spatial directions for arbitrary directions as described by Coscia et al. (2011). Here, we applied a weight of 0.1, which means that the penalty for purely vertical model gradients is 10 times lower compared with horizontal gradients, reflecting our general knowledge of the model (e.g., known by an ERT inversion result).

Figure 6 shows the resulting data misfits for two selected loops and all resistivity and inversion schemes. Loop rx 4 is located over the right slope of the buried valley structure and thus is anticipated to be affected most by the 2D resistivity structure, whereas rx 7 is over quasi-1D resistivity (Figure 5a). The misfits of the other measurements show similar features and a general agreement with those of rx 4 and rx 7.

In the given cases, we used a regularization parameter of $\lambda =250$. Higher values start to show systematic misfits, thus increasing the $\chi 2$ beyond 1. Except for the complex inversion of the 1D resistivity, all inversion cases fit the data within error bounds; that is, a $\chi 2$ of around 1.05 is achieved and the remaining differences between the observed data and the predicted data are randomly distributed. However, there are significant differences for the complex inversion using the 1D resistivity kernel (a final $\chi 2$ of 3.3), where significant systematic misfits are observed. For the layered-case complex inversion, we were not able to find a $\lambda $ that leads to a significantly better misfit; thus, we decided to present the results with identical regularization as for the other cases.

The corresponding models for water content and relaxation times are shown in Figures 7 and 8. The white lines correspond to the boundaries of the synthetic block model. However, they are only shown for visualization and are not used to construct the inversion mesh. All cases, apart from complex 1D case, are able to resolve the water content of the first aquifer in terms of layer boundaries to the aquitard as well as the lower edge and the asymmetric valley structure. The second aquifer below the valley/aquitard is less resolved. In particular, the boundary between the aquitard and second aquifer is difficult to detect. In terms of absolute values, the inversions are able to find the correct water content for the bottom layer. In general, the inversion of complex data shows improved resolution compared with those using rotated amplitude data.

Concerning the relaxation times (Figure 8), one sharp boundary associated with the first aquifer base is well-resolved. The absolute values of the two aquifers are close to the synthetic model. Only the clay layer cannot be detected because of the short relaxation time, as expected. The differences between the inversion results are negligible, with the exception of the complex 1D case.

As for the complex inversion of the 1D case, the resulting model shows strong artifacts that are not correlated to any of the expected structures together with unreasonably high water contents of up to 50% and low relaxation times of ≤40 ms. The artifacts are not surprising, considering the data misfit. Braun and Yaramanci (2011) already show that there are cases in which 1D approximations of the resistivity lead to erroneous kernels.

The effect on the inversion can be also seen when comparing different measurements (tx-/rx combinations). For example, the misfit for rx 4 (midpoint around $x=0\u2009\u2009m$, directly above the 2D structure) is affected more than rx 7 (midpoint around $x=75\u2009\u2009m$), which is half a loop diameter away from the 2D structure.

It seems that this is mainly a problem of the incorrectly estimated phase characteristics of the complex kernel function because part of the phase is defined through the resistivity structure. For the 1D case, the resistivity distribution introduces a phase in the kernel function that is significantly different from the kernel for the true blocky synthetic data set. The inversion tries to fit the data including the phase with changes in water content and relaxation times only, without success. This becomes clear when comparing the results and misfits to the rotated amplitude inversion. Because the amplitude approach circumvents the phase problem of the wrong resistivity estimation through the rotation of the data, the inversion is more stable and achieves much better results and misfits. However, the asymmetry in the second aquifer remains, even if the absolute values are more reasonable.

All inversion results show a general agreement when it comes to the relaxation times model. The $T2*$ time of the aquitard and the unsaturated zone cannot be detected because there is negligible signal from these units due to the short decay time or low water content, respectively. Apart from that, the two aquifers are clearly resolved in shape and absolute values. In general, the complex and amplitude inversions are in closer agreement than the water content models (Table 2). Even the complex inversion of the layered case shows a more reasonable model, compared to the water content. In general, the resistivity distribution seems to have a much lower impact on the relaxation times because the kernel directly targets the water content.

## CONCLUSION

We presented the open-source code COMET, which is able to model and invert surface NMR data. The model parameters, including resistivity, were allowed to vary in one dimension or two dimensions, and the loop configuration can be of arbitrary shape. Surface NMR data were inverted using rotated amplitudes or the complex data set using smooth (one dimension or two dimensions) or block inversion (one dimension). In addition, it is possible to use the forward operator for other inversion schemes such as Monte Carlo methods. We successfully used smooth resistivity distributions derived from ERT for the MRT inversion. The estimated models are very close to those using the original synthetic block resistivity distribution. Hence, the commonly used approximation of resistivity information as 1D or 2D blocks to set up an MRT inversion can be replaced by any arbitrary resistivity distribution. This is especially valuable for field cases. A comparison of rotated amplitude and complex inversions can hint to incorrectly estimated resistivities because a complex inversion is sensitive to the NMR phase. In cases in which the resistivity approximation is incorrect, the rotated amplitude inversion can achieve a good data fit; however, the water content model is prone to artifacts.

The COMET software package can handle (processed) MRS/MRT data in one dimension or two dimensions without the requirement of proprietary software. COMET enables control over all calculation steps and gives the user a powerful tool for sophisticated investigations, particularly for 2D resistivity and/or 2D NMR parameter distributions.

A structural coupling of NMR parameters and the resistivity could further enhance the image resolution of ERT and MRT. The total field approach from custEM could allow COMET to account for more complicated geometries such as topography, full-space applications, or arbitrarily shaped loops in three dimensions, for example, vertical loops at the wall of a tunnel. We are confident that the presented code can be picked up by the scientific community to expand MRS and MRT applications and simplify their documentation and reproducibility.

## ACKNOWLEDGMENTS

This research was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant MU 3318/3-1. We thank T. P. Irons and three anonymous reviewers as well as the associate editor F. Broggini for their extensive testing and care during the review process.

## DATA AND MATERIALS AVAILABILITY

The code is licensed under the copy-left GNU General Public License version 3 (GPLv3) and is freely available via GitLab (https://gitlab.com/Skibbe/comet) along with issue tracking and installation instructions. Detailed documentation is available on https://comet-project.readthedocs.io and includes the provided examples, tutorials to learn the main classes and their features, and a full API documentation.

### CODE FOR 1D MODELING AND INVERSION

Because most surface NMR applications consist of 1D soundings, we show an example of how to use COMET for a synthetic experiment with 1D layers (listing 1). This usually starts with the definition of the subsurface, in this case simple arrays for layer thickness, resistivity, water content, and relaxation times. This is used as input for instances of COMET (Figure 1).

from comet import pyhed as ph

from comet import snmr

import matplotlib.pyplot as plt

# info level for logger = console output

ph.log.setLevel(20)

# part 1: synthetic model (3-layer earth)

# thickness [m]

syn_thk = [10, 20]

# resistivity [Ohm*m]

syn_res = [1000, 200, 2000]

# water content [1]

syn_wc = [0.15, 0.25, 0.1]

# relaxation times [s]

syn_t2 = [0.2, 0.3, 0.1]

# part 2: loop configuration and resistivity

cfg = ph.config(rho=syn_res, d=syn_thk)

# 40*40 m square loop around origin [0, 0]

# max dipole distance 2m

loop = ph.loop.buildRectangle(

[[–20, –20], [20, 20]],

max_length=2, config=cfg)

# part 3: survey with earth and sounding

site = snmr.survey. Survey()

site.addLoop(loop) # only one loop

# magnetic field strenght + dir

site.setEarth(mag=48e-6, incl=60, decl=0)

# by default tx=0, rx=0 coincident loop

site.createSounding()

# part 4: setup 1D manager for inversion

# block inversion (or ‘smooth’)

# complex inversion (or ‘rotatedAmplitudes’)

mrs = snmr.modelling.MRS(

survey=site,

fid=0,

mtype=‘block’,

dtype=‘complex’)

# simulation of mrs data for given block model

mrs.simulate(

[syn_thk, syn_wc, syn_t2],

err=20e-9, num_gates=50, num_cpu=4)

# save kernel for later loading 4 lines above

mrs.kernel.save(‘kernel’)

# part 5: inversion with a fixed regularization

mrs.invert(lam=100)

print(‘final chi square:’, mrs.inv.getChi2())

# part 6: visualize results

mrs.showDataAndFit()

# comparison to syn model

mrs.showResult(syn=[syn_thk, syn_wc, syn_t2])

plt.show()

Listing 1. Synthetic modeling and inversion using the *COMET* toolbox.

In this case, we gather the resistivity in a primary configuration (simple class for the 1D model configuration). We assume a three-layer case with two finite layers of 10 and 20 m thickness and an infinite layer. The middle layer bears more water (25% compared with 15% and 10%) and is less resistive ($200\u2009\u2009\Omega m$ compared with 1000 and $2000\u2009\u2009\Omega m$). Further, we define a 40 m square loop that is used as transmitter and receiver, which also holds the resistivity distribution.

A bare survey class instance is used as main object to gather the above-defined information. The survey class manages the different loops, soundings, and data sets as well as the earth’s magnetic field. It is also used to create sounding class instances representing FID experiments defined through given tx and rx indices.

The MRS class is a derivative of the modeling or manager class for a 1D magnetic resonance sounding with either block or smooth NMR parameter distributions. It gets a reference to a survey class as well as the index of the corresponding FID (in this case, the first and only number 0).

The simulation is started via a simple method call. Up to this point, we only initialized the different objects, but did not calculate anything. The MRS class tries to create a synthetic data set and asks the kernel class for a kernel matrix. The corresponding kernel class has no matrix yet and starts calculating one; however, it needs magnetic fields. Through the given tx index (implicit through createSounding(…)) it finally reaches the loop class, which then starts computing a magnetic field. Finally, the inversion is carried out with respect to the given parameters and desired model and data type by calling invert(…).

This whole principle of supply-on-demand is used throughout the code. However, it is also possible to calculate the magnetic field kernel first, save it, and load it at any later stage so that the kernel computation is not carried out if only an inversion parameter is changed.

Note that the simulated data are stored in the FID instances, hence the survey class. For saving the data, the survey class method save(…) can be called. The MRS instance only provides methods for modeling and inversion; it does not contain any information other than inversion parameters. Moreover, it also provides several convenience functions for visualization based partly on the general-purpose plotting library matplotlib and partly on pyGIMLi.

### CODE FOR KERNEL CALCULATION IN TWO DIMENSIONS

This part describes the working steps for the 2D example and is divided into the individual steps visible in Figure 5. First, we import the different packages in our script (Listing 2). This enables us to use their functionalities in the other code snippets. Additionally, we do some definitions of the earth’s magnetic field and the background resistivity. Finally, we import a 2D parameter mesh with a suited resistivity vector, for example, that saved after an ERT inversion.

# import main parts of COMET

from comet import pyhed as ph

from comet import snmr

# import numpy and pygimli libraries

import numpy as np

import pygimli as pg

# earth magnetic field

earth = snmr.survey. Earth(

decl=0, incl=60, mag=48000e-9)

# loop configuration: 1D earth + frequency

# e.g. halfspace with 1000 Ohmm

loop_config = ph.config(

rho=[1000.],

f=earth.larmor)

# import 2D parameter mesh from ERT inversion

paramesh2d = pg. Mesh(‘invmesh.bms’)

# import associated resistivity vector

res = np.load(‘anom_res.npy’)

Listing 2. Import of used modules and definition of magnetic field and resistivity.

Next, we create an FEM mesh suited for a given ETRA layout (Listing 3). This snippet represents the Python code used to get from Figure 5b and 5c. The mesh and resistivity distribution are fed into the temporarily created loop class (mesh_loop). Note that a boundary is appended to ensure boundary conditions. The creation of a suitable FEM mesh as well as the interpolation of the resistivity from two dimensions to three dimensions is mostly done automatically. The user can have a look at the resulting meshes, for example, using ParaView (Ahrens et al., 2005).

# etra survey, small edge length = 50 meter

# using 8 loops in total

# configuration with res and frequency is given

loops = ph.loop.buildEtraSurvey(

50, num_loops=8, config=loop_config)

# simple merge measurement geometry

# of all loops

mesh_loop = ph.loop.mergeLoops(loops)

# set 2D mesh and geometry

mesh_loop.setParaMesh2D(

paramesh2d, append_boundary=True,

anomaly=1./res) # need conductivity here

source_poly = ph.loop.buildEtraPoly(

-100, 100, 50)

# finally create FEM mesh and export

# additional VTK for quality control

# define refinement box

mesh_loop.createFEMMesh(

savename=‘fem_mesh’,

source_poly=source_poly,

exportH5=True, box_x= [–125, 125],

box_y= [–125, 125], box_z=–50,

inner_area_cell_size=2,

subsurface_cell_size=35.)

Listing 3. Building a custom FEM mesh for complex loop configurations.

The calculation of the magnetic field (Figure 5c and 5d) is done using the previously generated FEM mesh (listing 4). To this end, we associate the configuration and the resistivity distribution for every loop, before we compute the primary and secondary fields that are added and saved.

# loop over all loops to be calculated

for lp, loop in enumerate(loops):

# set 2D ERT inversion mesh with boundary

loop.setParaMesh2D(

‘invmesh.bms’,

anomaly=1./res,

append_boundary=True)

# import FEM mesh

loop.setFEMMesh(‘fem_mesh.bms’)

# prepare primary field for secondary

# field calculation on 48 cores

loop.prepareSecondaryFieldCalculation(

num_cpu=48)

# calculate secondary field

loop.calculateSecField(num_cpu=48)

Listing 4. Calculating magnetic field for the used loops.

Finally, the 2D kernel function is calculated on the 2D ERT parameter mesh (listing 5). Previously cached loop class instances (plus magnetic fields) and meshes are imported for this purpose.

All of the main classes (survey, loop, FID, and kernel) can be saved to disk and loaded at a later stage.

# We want to gather everything in this bare

# survey class. Later we can simply import the

# whole thing

site = snmr.survey. Survey()

# looping over the same loops

for lp, loop in enumerate(loops):

# add loop to survey class

site.addLoop(loop)

# create FID sounding with

# first loop as tx

# eight loops lead to eight

# different FID’s = measurements

site.createSounding(tx=0, rx=lp)

# load loopmesh (cached internally

# to save time)

site.loadLoopMesh(‘femmesh.bms’)

# give to survey

site.setEarth(earth)

# do for all FID’s

for i in range(len(site.fids)):

# generate kernel class for a 2D

# kernel for fid index i

kern = site.createKernel(fid=i, dimension=2)

# define kernel base and refinement

kern.set2DKernelMesh(‘invmesh.bms’, order=2)

# calculate kernel in 2D slices in parallel

kern.calculate(num_cpu=48)

kern.save(‘kern_{}’.format(i))

Listing 5. Calculating the kernel for the used ETRA setup.