# An overview of a highly versatile forward and stable inverse algorithm for airborne, ground-based and borehole electromagnetic and electric data

Esben Auken, Anders Vest Christiansen, Casper Kirkegaard, Gianluca

Fiandaca, Cyril Schamper, Ahmad Ali Behroozmand, Andrew Binley, Emil Nielsen, Flemming Effersø, Niels Bøie Christensen, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Esben Auken, Anders Vest Christiansen, Casper Kirkegaard, Gianluca Fiandaca, Cyril Schamper, et al.. An overview of a highly versatile forward and stable inverse algorithm for airborne, ground-based and borehole electromagnetic and electric data. Exploration Geophysics, 2018, 46 (3), pp.223-235. 10.1071/EG13097. hal-03867420

HAL Id: hal-03867420

https://hal.science/hal-03867420v1

Submitted on 19 Jun 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

An overview of a highly versatile forward and stable
inverse algorithm for airborne, ground-based and borehole
electromagnetic and electric data

Esben Auken*¹, Anders Vest Christiansen¹, Casper Kirkegaard¹, Gianluca

Fiandaca¹, Cyril Schamper², Ahmad Ali Behroozmand¹, ⁷, Andrew Binley³, Emil

Nielsen⁴, Flemming Effersø⁵, Niels Bøie Christensen¹, Kurt Sørensen¹,⁵, Nikolaj

Foged¹ and Giulio Vignoli¹, ⁶

*Corresponding Author: Department of Earth Sciences, Aarhus University, C.F.

Møllers Allé 4, DK-8000 Aarhus C, Denmark

¹ Aarhus University, Denmark

² Sorbonne Universités, UPMC Univ Paris 06, France, formerly at 1.

³ Lancaster University, GB

⁴ Danish Technical University, Denmark

⁵ SkyTEM Surveys, Denmark
21 ⁶ Geological Survey of Denmark and Greenland (GEUS), Denmark, formerly at

22 1

23 ⁷ Stanford University, USA, formerly at 1

24
25

26

ABSTRACT

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

Geophysical methods for high resolution studies of the near surface have become a powerful tool for many different modern applications of widely varying scope and scale. Often geologists, engineers or hydrologists use a custom combination of different geophysical techniques for their particular purpose involving dedicated inversion algorithms. This is inconvenient and makes direct result comparison difficult due to different model formulations, approximations, regularization, lack of uncertainty estimates etc. We present an overview of a mature, robust and general algorithm providing a single framework for the inversion of most electromagnetic and electrical data types and instrument geometries. The implementation mainly uses a 1D earth formulation for electromagnetics and magnetic resonance sounding responses, while the geoelectric responses are both 1D and 2D and the sheets response models a 3D conductive sheet in a conductive host with an overburden of varying thickness and resistivity. In all cases focus is put on delivering full system forward modeling across all supported types of data. Our implementation is modular, meaning that the bulk of the algorithm is independent of data type, making it easy to add support for new types of data. Having implemented forward response routines and file I/O for a given data type provides access to a robust and general inversion engine. This engine includes support for mixed data types, arbitrary model parameter constraints, integration of prior information and calculation of both model parameter sensitivity analysis and depth of investigation. We present a review of our implementation and methodology and show four different examples illustrating the versatility of the
algorithm. The first example is a Laterally Constrained joint Inversion (LCI) of
surface time domain induced polarization (TDIP) data and borehole TDIP data,
the second example shows a spatially constrained inversion (SCI) of airborne
transient electromagnetic (AEM) data, the third example is an inversion and
sensitivity analysis of Magnetic Resonance Sounding (MRS) data, where the
electrical structure is constrained with AEM data. The fourth example is an
inversion of AEM data, where the model is described by a 3D sheet in a layered
conductive host.

INTRODUCTION

A wide range of near surface electric and electromagnetic geophysical
measurement techniques are routinely employed for vastly different purposes
within a variety of disciplines in modern science and engineering. Combining
advanced forward modeling with sophisticated inversion schemes allows for
obtaining accurate information on electric resistivity properties of geological
layers, however, such algorithms are typically limited to the modeling of a single
or a few types of data. Examples of this include

• 1D inversion of frequency domain ground conductivity meter data
(Santos 2004),

68 • airborne frequency domain data in 1D (Sengpiel and Siemon 2000),

69 • versatile inversion of frequency domain EM data in 1D (Oldenburg and Jones 2011b),

71 • 1D holistic inversion of airborne frequency- (Brodie and Sambridge 2006a) and time domain EM data (Brodie and Sambridge 2006b),
72
73

* 1D inversion of frequency and time domain data (Christensen and Auken 1992),

74

75

* 1D approximate inversion of time domain data (Tartaras et al. 2000; Christensen 2002; Macnae et al. 1991),

76

77

* 1D inversion of multiple types of time domain data (Oldenburg and Jones 2011c),

78

79

* 3D inversion of airborne time- and frequency domain data (Newman and Alumbaugh 1997; Cox, Wilson and Zhdanov 2010, Oldenburg, Haber and Shekhtman 2013),

80

81

82

* 2D inversion of resistivity (ERT) and induced polarization data (IP) (Loke and Barker 1996; Loke, Chambers and Ogilvy 2006; Oldenburg and Li 1994, Fiandaca et al. 2013b),

83

84

85

* 3D inversion of ERT data (Günther, Rücker and Spitzer 2006; Rücker, Günther and Spitzer 2006),

86

87

* 3D inversion of ERT data with IP (Yoshioka and Zhdanov 2005; Loke 2011; Oldenburg and Jones 2011a),

88

89

* 3D inversion of magnetotellurics data (Gribenko, Wan and Zhdanov 2010),

90

91

* inversion of Magnetic Resonance Soundings (MRS) recovering water content (Mohnke and Yaramanci 2002; Mueller-Petke and Yaramanci 2010; Hertrich, Braun and Yaramanci 2005)

92

93

94

Implementing just the base of such inversion algorithms can be a tremendous task in itself and often resources will limit the stage of development that can be reached. For this paper we present a more unusual research algorithm, AarhusInv, in the sense that it supports a broad spectrum of data

95

96

97
types. AarhusInv has further been allowed to reach a point of production maturity, while at the same time being very actively developed for various research projects. Reaching this stage has to a large degree been made possible by support from the ongoing Danish national mapping of groundwater resources (Thomsen, Søndergaard and Sørensen 2004; Møller et al. 2009). This project was initiated in the late 1990's and has provided stable funding and extensive production use of the algorithm, including inversion of ERT data, PACES (Sørensen et al. 2005), groundbased TEM and airborne SkyTEM data (Sørensen and Auken 2004). Furthermore, there have been many contributions on various parts of the algorithm, hence the long list of co-authors on this paper.

The main distinction to other algorithms reported in the literature is how different high accuracy forward models, spanning multiple data types, are brought together in a common inversion framework. This framework is not only optimized for production but also flexible for research, by supporting arbitrary spatial constraints and full integration of a priori information on any model parameter.

AarhusInv is freely available for university research institutions. Over time it has been used for the inversion of an extensive amount of vastly different datasets, collected over all continents of the world. Some of the most prominent include an 18,000 line km VTEM survey conducted over the Okavango delta, Botswana (Podgorski et al. 2013), a 34,000 line km SkyTEM survey in the Broken Hill Area of Australia in 2009, 14,000 line km in India in 2013 and a 1,000 line km TEMPEST survey on Eyre peninsula of South Australia (Auken et al. 2009a). For the Danish groundwater mapping project it has been used for inversion of approximately 40,000 groundbased TEM soundings, 25,000 line
123 kilometers of SkyTEM and several thousand ERT profiles. We estimate that the
124 algorithm has inverted more than 400,000 line kilometers of airborne data since
125 2005, including data from the Resolve system (Fugro Inc.), DigHEM (Fugro
126 Inc.), SkyTEM (Sørensen and Auken 2004), VTEM (Witherly, Irvine and
127 Morrison 2004), AeroTEM (Balch, Boyko and Paterson 2003) and TEMPEST
128 (Lane et al. 2000).

129 In this paper we provide a review of the numerical implementation and
130 demonstrate its flexibility and capabilities through examples. A core component
131 is a constrained model framework which has proven successful for many
132 different geologies, ranging from mapping of paleo channels in sedimentary
133 environments to mapping of perched aquifers in a weathered volcanic geology
134 (Jørgensen and Sandersen 2006; d'Ozouville et al. 2008; Danielsen et al.
135 2003; Reid, Munday and Fitzpatrick 2007; Siemon et al. 2002). We begin by
136 describing the design of the algorithm along with a list of supported data types
137 and instrument geometries. This is followed by a review of the forward and
138 inverse modeling schemes. Finally four field data examples illustrate the
139 algorithms versatility to operate on data collected on the ground, in the ground
140 and in the air.

# METHODOLOGY

142 **Modular algorithm layout**

143 Since the origin of AarhusInv in 1995 - 1997 its feature set has evolved
144 significantly, but in its simplest form the functionality is still to invert a set of
145 geophysical soundings for a set of layered 1D models connected through
146 constraints. The algorithm was formerly known as em1dinv, but with the
147 introduction of multidimensional responses we changed the name to AarhusInv
148 in 2012. The code base is in Fortran 95/77 and the program itself is a command
149 line application with no user interface and with all input/output conducted
150 through ASCII files. This structure is flexible and very well suited for functioning
151 as an integrated inversion engine, as in the case of the front end GUI
152 applications Aarhus Workbench (Auken et al. 2009c), EMMA (Auken et al.
153 2002) and SPIA (www.hgg.au.dk). The algorithm is implemented in a general
154 modular manner as conceptually outlined in Figure 1. This figure is
155 supplemented by Table 1, providing full details on supported data types,
156 measurement geometries, and key references.

157 Starting from the top in Figure 1 a general model input file of flexible
158 format is first read. This model file contains a starting model definition, specifies
159 arbitrary regularization constraints and also holds prior information. Following
160 model file input is the reading of data files in the data file input module. These
161 files are encoded in data type specific formats, which are transparently handled
162 by the modules internal branching to the relevant sub modules. Note that mixing
163 of any number of different data types is supported, allowing for unrestricted joint
164 inversions. Having filled the internal data structures from input files these
165 structures are passed on to the general inversion module which starts the
166 iterative inversion process. During the iterative course of minimizing the
167 residual the inversion module calls the forward response module (see Appendix
168 A) for the calculation of forward response updates and derivatives. The forward
169 response module also handles data type specific branching transparently,
170 simply performing the requested calculations and returning the result to the
171 inversion module. Internally, the forward response module makes data type
172 dependent calls to a library of common routines. Having iteratively solved the
173 inverse problem (see next paragraph and Appendix B) the result is sent to the
174 general model analysis module, which calculates a model parameter sensitivity
175 analysis. Finally the depth of investigation (DOI) is calculated based on the
176 actual model output from the inversion process and a recalculated sensitivity
177 (Jacobian) matrix (Christiansen and Auken 2012). Following the DOI calculation
178 a result file is finally written in a general file format.

179 Forward modeling

180 From a forward modeling point of view AarhusInv provides a wide range
181 of options for simulating a given measurement system. This is particularly
182 important for TEM data, as neglecting the influence of system parameters such
183 as geometry, waveform or filters can lead to serious modeling error as
184 described for TEM by (Christiansen, Auken and Viezzoli and Fiandaca et al.
185 (2012) for DC/IP. The extent to which the algorithm allows for full system
186 forward modeling is relatively unique and makes it possible to compare model
187 results from different data types as directly as possible. Key references for the
188 forward modeling are given in the caption of Table 1. It is outside the scope of
189 this paper to describe all these different solutions and therefore we have in
190 Appendix A chosen to give descriptive details on the 1D EM solutions which are
191 also the most used parts of AarhusInv.

192 The forward modeling is an impended parallel using OpenMP and scales
193 close to linearly running on 1 – 64 processors. All vectors and matrices are
194 furthermore implemented sparsely making the memory consumption relatively
195 low even for very large inverse problems. It is out of the scope of this paper to
196 describe the details of this implementation and we therefore refer to Kirkegaard
197 and Auken (2014) for all details.

198

199 **Inverse modeling**

200 The mathematical formalism of the inverse modeling scheme follows the
201 established practice by Menke (1989a). Details are provided by Auken et al.
202 (2005) and Appendix B gives in-depth details to the implementation of the
203 inversion. Here, we focus on presenting the flexibility and capabilities of the
204 algorithm. Essentially, we use a Gauss-Newton style minimization scheme with
205 a Marquart modification (Marquart 1963) to find the set of model parameters
206 that minimize the L₂ misfit with respect to observed data, regularizing and prior
207 information. All data sets are inverted simultaneously, minimizing a common
208 object function. Consequently, the output models are balanced between the
209 constraints, the physics and the data. Model parameters that have little
210 influence on the data will be controlled by the constraints, and vice versa. The
211 use of, for example, lateral or spatial constraints on a collection of models
212 allows for information to propagate across a survey area.

213 Model constraints are applied between any two parameters of the same
214 type by specifying the variance of the difference between the two and the
215 uncertainty of prior information is similarly specified by the variance of the given
216 prior value. The AarhusInv implementation itself is unbiased in terms of how
217 constraints should be defined and thus leaves this to be specified by the user
218 in the input model file. This format allows for arbitrary constraints linking any
219 model parameters, i.e. a set of 1D models can be linked by 3D constraints in
220 any way the user sees fit. Typically, we apply constraints in the form of 1D
221 laterally constrained inversion (LCI, Auken et al. 2005) and 1D spatially
222 constrained inversion (SCI, Viezzoli et al. 2008) for producing quasi 2D and
223 quasi 3D models, respectively. The constraining formulation is flexible and
224 allows for constraints and prior information on both primary and secondary
225 model parameters. We define primary model parameters as the actual
226 parameters of the layered model, e.g. resistivities and thicknesses, whereas
227 secondary model parameters are any parameters that can be derived from
228 linear combinations of the primary parameters, e.g. depths and elevations
229 (Auken et al. 2005). This can be utilized for example when integrating prior
230 information from a borehole, where the known parameter is typically the depth
231 of a layer interface and not a layer thickness.

232 For the inversion of airborne data we allow including the instrument
233 altitude as a model parameter, as will be discussed in more detail in the airborne
234 example. In the case of helicopter data the pilot attempts to follow the terrain
235 topography and maintains a relatively constant altitude, whereas in fixed wing
236 surveys the aircraft typically operate at almost constant elevation. In this case
237 it is desirable to supply lateral constraints linking the modeled instrument
238 elevation of neighboring models, instead of linking the primary model
239 parameters in the form of height above terrain. Table 2 shows which additional
240 model parameters are included in the inversion for certain airborne systems.
241 These parameters are all in excess of the primary and secondary parameters
242 of resistivity, thickness and depth.

243 In conjunction with any inversion result a model parameter analysis can
244 be calculated, obtaining estimates of the resulting model uncertainty. This is
245 done by calculating the covariance matrix of the resulting model of the inversion
246 (Tarantola and Valette 1982a, Auken and Christiansen 2004), which provides
247 a linearized error estimate. Since we perform our inversion in logarithmic model
248 space for numerical stability, the values obtained from such a sensitivity
249 analysis can conveniently be regarded relative measures of uncertainty.
250 Absolute error estimates from logarithmic space are translated into a standard
251 deviation factor in linear space, such that 1.0 is equal to perfect resolution and
252 1.1 corresponds to a standard deviation of approximately 10%. We find that
253 when such linearized error estimates become much larger than unity they are
254 best viewed as guidelines, however. The use of such error estimates will be
255 illustrated in the following examples.

256 We finalize the inverse modeling procedure by calculating the DOI for
257 the resulting output models, an operation which relies on a reweighting of the
258 Jacobian matrix. The computations employ a global and absolute sensitivity
259 threshold value, which has been tuned for operating in logarithmic model/data
260 space (Christiansen and Auken 2012). For a given model, the DOI calculations
261 consider only the parts of the Jacobian related to observed data, implying that
262 the effect of lateral, spatial or vertical model parameter constraints and a priori
263 information is not included. With this type of DOI estimate it is possible to judge
264 when the information in a model is driven by data or heavily dependent on the
265 starting model or specifics of the regularization.

266 Limitations due to dimensionality of the forward response and inversion

267 The functionalities described so far illustrate the general capabilities of
268 the algorithm at its current stage of development, but obviously its inherent
269 limitations should also be considered.
First of all the algorithm base began as a 1D implementation and the whole suite of EM methods are still limited to 1D forward modeling with the exception of 3D sheets. For DC and IP the forward algorithms are present in both 1D and 2D versions and a 3D implementation is close to being complete. Solving the physical problem in 1D can be a limitation for applications such as mineral exploration, whereas studies of sedimentary settings fits well within a 1D modeling envelope (eg. Auken et al. 2008, Viezzoli et al. 2010).

Secondly, while being implemented as general as possible, our regularization and inversion schemes are fundamentally based on a static formulation. This means that it is not possible to experiment with e.g. Tikhonov-style variable regularization parameters. However, the least squares (L2-norm) solution can be changed to use the L1-norm or the sharp boundary formalism by Zhdanov, Fang and Hursán (2000) and Vignoli et al. (2013).

## EXAMPLES

In the following we show four examples demonstrating the capabilities and versatility of the AarhusInv algorithm. In the examples we show inversions of data collected on the ground, in the ground and in the air. We also illustrate the use of spatial-, lateral- and vertical constraints. The examples are chosen to illustrate widely varying scales of interest, ranging small scale engineering type of problems to large regional scale surveys. We further illustrate the use of prior information and provide a synthetic example of joint inversion of mixed type data. Finally we provide a synthetic example inverting an airborne data set with a 3D sheet model.
293 **On and in the ground: DC and time domain induced polarization**

294 For our first example we present a novel application of joint 1D LCI
295 inversion of surface time domain induced polarization (TDIP) data and borehole
296 TDIP data. The TDIP method is a natural and efficient extension of standard
297 DC multi electrode profiling by simply adding the measurement of the time
298 dependence of discharge after an injection current is turned off. Modeling the
299 complete time decay of each data point in the Cole-Cole formulation (Pelton et
300 al. 1978) allows for extracting additional independent model parameters, often
301 making it possible to discriminate earth structures with an otherwise identical
302 signature. Apart from resistivity, each layer of an IP model includes the
303 parameters τ, C and M₀. Here m₀ is the magnitude of the chargeability when the
304 injection current is turned off at t=0, τ is the decay time constant and C is a
305 constant controlling the frequency dependence of the response. The TDIP
306 implementation is described in full detail by Fiandaca et al. (2012) and includes
307 modeling of the full injection current waveform and low pass filters, which can
308 be a cause of serious modeling error when neglected. Compared to the
309 methodology presented by Fiandaca et al. (2012) we add in this example
310 borehole TDIP data for a joint inversion of two data types sharing a common
311 model space, in order to demonstrate the versatility of the algorithm.

312 The investigated area is a former Danish landfill active from the 1950's
313 to the 1980's. During this period the landfill has been almost complete
314 unregulated and it was further established without any kind of membranes,
315 leachate capture or isolation systems. Previous studies show that oils and
316 chemical waste from the landfill is contaminating the area by water seepage
317 through the landfill, but the extent of the pollution is largely unknown. Additional
318 details of the survey area are provided by Gazoty et al. (2012). Figure 2 show
319 the dataset consisting of DC/IP data collected using a gradient array protocol
320 and an electrode spacing of 5 meters, as well as DC/IP data from a pole-pole
321 El-log (Sørensen and Larsen 1999) collected in a borehole around the center
322 of the profile. For the logging a measurement was made for every half meter
323 down to 24 meters. This combined dataset was 1D LCI inverted for models of
324 21 layers, discretized in fixed boundaries and utilizing lateral and vertical
325 smoothness constraints.

326 On the model sections in Figure 2 the position of the borehole is
327 indicated near the center. In the right part of the sections further results for 2
328 drillings are superimposed for comparison. From these boreholes it is tempting
329 to conclude that the waste and pollution is likely to be localized to the resistive
330 structure in the vicinity of the boreholes. However, when including the
331 chargeability (m₀) section in the interpretation it becomes clear that the extent
332 of the waste body can possibly be much greater. In fact the results of the full
333 investigation of 13 profiles and 15 boreholes show a very strong correlation
334 between buried waste and chargeability, which is undetectable in terms of
335 resistivity alone. Including log data in a joint inversion contributes
336 complementary information due its perpendicular plane of sampling, introducing
337 a completely different frequency domain kernel since the electrodes are placed
338 in the ground. The result of including this complementary and very detailed
339 vertical information is for all model parameters to become well determined at
340 this position as seen in the right part of Figure 2. In turn the entire profile benefits
341 from this added information, since the use of LCI constraints enable lateral
342 migration of information.
343 **In the air: Transient electromagnetic**

344 The vast majority of data being inverted by AarhusInv is airborne data from
345 helicopter and fixed wing systems. To illustrate this key application we show
346 result of an SCI inversion of SkyTEM data (Sørensen and Auken 2004). The
347 SkyTEM system uses dual transmitter moments to provide unbiased time
gates from as early as 6 µs up to 10 ms (Schamper, Auken and Sørensen
348 2014). This allows for accurate forward modeling, including the effects of the
349 waveforms of the dual transmitter moments, low pass filters, front gate and
350 finite width gates. For the inversion we use measurements of the vertical
351 component of the secondary field only and include instrument altitude as a
352 model parameter. In the more general case of fixed wing type geometry the
353 algorithm allows for inversion of all field components, including receiver pitch
354 and roll as model parameters (Auken et al. 2009b). We have also
355 implemented inversion for receiver coil response model parameter, allowing
356 for the use of data from time gates as early as 6 µs (Schamper, Auken and
357 Sørensen 2013)
358

359 Our example dataset of 3250 line km was collected in 2009 covering an
360 area of approximately 820 km² intersected by the Danish/German border,
361 having virtually no topography and covering the coast to the North Sea. The
362 survey was conducted to provide a better overall understanding of the
363 hydrogeological setting in the area and to map the fresh/salt water interface in
364 the costal region. Following data acquisition, all couplings due to man-made
365 installations were removed from the raw data set and further processed using
366 the methodology of Auken et al. (2009c). The result is a data set of stacked
soundings for each transmitter moment approximately every 25 meters. This data set was jointly inverted for both moments using 19 layer 1D models with fixed layer boundaries and constrained using the SCI methodology. We show the results of the inversion in Figure 3; sub figure (a) shows the flight lines of the dataset, (b) the depth of investigation, (c) the modeled resistivity at a depth of 74-88 meters and in (d) the model parameter analysis of (c). From resistivity maps (c) it is possible to identify important features of both regional- and more local scale if zooming into the details. This is illustrated by marking the extent of two buried valley structures and a large area of salt water intrusion. The geological interpretation of the dataset including cross sections and average resistivity maps is discussed in detail by Jørgensen et al. (2012).

## 379 **Mixed data types: Magnetic Resonance Sounding and Airborne TEM**

For this last example we show a feasibility study involving forward modeling, inversion and sensitivity analysis of Magnetic Resonance Sounding (MRS) data in combination with airborne TEM data. The example is intended to illustrate the capabilities of the algorithm on mixed type datasets, utilizing multiple modules of the algorithm simultaneously.

The synthetic model of our study is outlined in Table 3 and represents a 20 meter thick aquifer situated at a depth of 50 meters, covered by dry sands, a till layer and defined at depth by a thick layer of clay. For this synthetic model we consider two types of synthetic data: an airborne TEM sounding and a ground-based MRS sounding. The TEM sounding simulates a SkyTEM system utilizing a 300 m² transmitter loop situated at an altitude of 30 meters, acquiring data for time gates ranging from 10 µs to 1 ms. For the inversion the simulated
392 data is assigned a uniform noise level of 3%, but no actual perturbation is
393 performed. In case of the MRS part of the dataset we consider a simulation
394 utilizing a 100 m x 100 m square loop. We further set the magnetic field of the
395 earth to 49,300 nT at an inclination of 70 degrees and a declination of +2
396 degrees and utilize pulse moment values ranging from 0.11 to 15.0 As. The
397 synthetically generated sounding was assigned and perturbed for a base noise
398 level of 20 nV, with an additional uniform noise contribution of 3%.

399 Traditionally MRS data has been inverted in what can be referred to as
400 step-wise inversion schemes (e.g., Günther et al. 2011; Mueller-Petke and
401 Yaramanci 2010; Behroozmand et al. 2012b). This type of inversion starts by
402 inverting a TEM or DC sounding for a resistivity model, which is then assumed
403 to be the true model in a subsequent inversion for MRS model parameters
404 (water content ratio W and relaxation time T2*). In other words the MRS kernel
405 function, which is a function of the resistivity model, is assumed fixed.
406 Behroozmand et al. (2012a) proposes a joint inversion methodology where the
407 resistivity model parameters are free in the inversion, implying that the MRS
408 kernel has to be updated for each iteration of the inversion requiring a fast
409 numerical implementation of the MRS kernel. In the following we compare the
410 resolution capabilities of these two methods.

411 In Figure 4 the model results of a traditional step-wise inversion of the
412 simulated data is shown. The figure features dashed blue lines for indicating
413 68% confidence intervals, as obtained from a sensitivity analysis. It is clear that
414 the upper- and lower-most boundaries of the resistivity model are very well
415 resolved, whereas the layering in between these boundaries is completely
416 undetermined. If a noise perturbation had been applied to the TEM data, the
417 obtained model result would have been very far from the true model for the two
418 highly resistive layers. Turning to the MRS model we see that the water content
419 of the synthetic model is well resolved for the top two layers, but underestimated
420 for the lower lying water bearing sand layer with the true value out of range of
421 the models 68% confidence intervals. This mismatch is caused by the slightly
422 wrong layer thicknesses fixed in the resistivity model, which in turn forces the
423 modeled water content to compensate for the MRS kernel function being
424 slightly off. For the bottom layer the inverted water content is found to be even
425 further from the true model, but in this case the deviation can to a larger degree
426 be attributed to a high degree of uncertainty as it is closer to the depth of
427 investigation. In general, T2* values are found to be reasonably well resolved
428 for all layers.

429 In Figure 5 the result of the joint inversion scheme is shown. The first
430 thing to note is how all layers of the resistivity model have become well
431 resolved, also in the case of the two resistive layers. Determination of layer
432 thicknesses has improved in particular, due to the combined, complementary
433 information provided by the SkyTEM and MRS signals. In the case of MRS
434 model parameters this point can be seen from the water content model being
435 very accurately reproduced for the top three layers, while the improvement for
436 the deepest lying layer is less pronounced as it is closer to the depth of
437 investigation. For the T2* model the impact of the coupled scheme is very little.

438
Airborne TEM and sheets inversion

In the AarhusInv code, the forward modeling of thin sheets is based on
the algorithm developed by Zhou (1989) with the mathematical formulation from
Weidelt (1983). This algorithm calculates the response from a number of
439
440
441
442 arbitrarily located thin sheets in a conductive layered host. The sheets have
443 both inductive and channeling current modes. If several sheets are used they
444 will be inductively coupled. To speed up the computation the algorithm uses
445 OpenMP to take advantage of multi-processor computer architecture, and an
446 efficient solver has been added to compute the scattered fields for inversion of
447 airborne datasets where many source-receiver locations need to be computed.
448 Furthermore, the algorithm has been extended to include both dipole and finite
449 loop sources, and it calculates responses both in the time and frequency
450 domain as stated in Table 1 and Table 2.

451 We will illustrate the use of this part of AarhusInv on synthetic transient
452 AEM data as shown in Figure 6. The model has a 20 m conductive overburden
453 of 50 Ωm, covering a resistive halfspace of 1000 Ωm. The sheet is in the
454 resistive host and has a conductance of 50 S, a top depth of 30 m, a maximum
455 depth of 175 m, and both its dip and strike angles are 45° (Figure 6a). The
456 simulated AEM dataset has a line spacing of 100 m with a sounding spacing of
457 20 m along the lines (Figure 6b) and z-component data only. The flight altitude
458 is 30 m and recorded times span from about 10 µs to a few ms depending on
459 the noise. Random noise at a slope of t⁻¹/² and a noise level of 1 nV/m² at 1 ms
460 was added to the data (Schamper et al. 2014). Figure 6c illustrates the fit
461 obtained after inversion with some of the recorded times along the central
462 profile. The starting model parameters are indicated in the figure. The
463 parameters in the inversion are sheet conductance, dip, strike, layered model
464 parameters and flight altitude. In the present case, the sheet is well determined
465 with the TEM method and the final estimated model is well recovered and
466 almost identical to the true model as is also indicated by the parameter determination (not shown here).
467

468 CONCLUSIONS

469 We have presented an overview of a proven, versatile and flexible
470 production inversion algorithm, capable of inverting most electric and
471 electromagnetic data types and supporting most instruments normally used
472 within the field of near surface geophysics. The algorithm is freely available for
473 academic use and provides the user with freedom to specify arbitrary 3D
474 regularization constraints, include any amount of prior information and invert for
475 mixed type data sets. Our implementation focuses on consistently accurate
476 system forward modeling independent of data type and instrument, which not
477 only ensures accuracy, but can also be regarded a prerequisite for joint
478 inversion of data collected using instruments of different transfer function. The
479 implementation is modular, meaning that the bulk of the algorithm is
480 independent of data type, making it very easy to add support for new types of
481 data. Calculation of both model parameter sensitivity analysis and depth of
482 investigation is a key feature and it is handled automatically regardless of
483 configuration and data types.
484

The versatility of the algorithm is illustrated by four different examples.
484
485 First, we showed the joint inversion of a resistivity/time-domain IP data set with
486 electrodes both located in the ground and on the surface. This example
487 provides an illustration of how including Cole-Cole parameters in the model
488 greatly improves the delineation of a former landfill. Second, we have shown
489 results of a large scale airborne TEM survey where the combined information
490 from maps of average resistivity, depth-of-investigation and model parameter
analysis improves the understanding of the underlying geology. Thirdly, we
have shown a synthetic example where we jointly invert Airborne TEM data with
ground-based MRS data. These two model-spaces are only vaguely related,
but we show how the mutual information significantly improves the resolution of
both the resistivity- and the water content model. In the fourth example we
showed a synthetic example inverting Airborne TEM data with a 3D sheets
model.

498

499

ACKNOWLEDGEMENTS

Over the years the development of the code has been supported by a
number of different research projects. The most prominent has been the HOBE
Center of Excellence (Villum Kahn Rasmussen Foundation) and the four
projects RiskPoint, NICA, CO2GS and HyGEM, all larger research projects
funded by the Danish Counsel for Strategic Research. Also the Nature Agency
under the Danish Ministry of the Environment has greatly supported our work
as an integrated part of the national groundwater mapping program. Lastly, we
have to thank numerous students who had to suffer long nights with new
implementations not always working as expected in the first beta-releases.

509
# APPENDIX A – 1D EM FORWARD CALCULATION

The type of problem we consider consist of a given instrument source, receiver and layered half space for which to solve Maxwell's equations. We do this in order to calculate some or all of the electromagnetic field components, or their time derivatives, at the position of the receiver. This type of calculation is conveniently performed in the frequency domain and can be formulated to exploit symmetry by choosing the gauge leading to the Schelkunoff potential formulation. In the following we closely follow Ward and Hohmann (1988). In the Schelkunoff formalism we operate on the electric and magnetic vector potentials **A** and **F** from which the actual electromagnetic fields can be determined by means of differentiation. This formulation is desirable since the vector potentials point in the direction of a source with a convenient geometry, allowing for calculation of all EM field components from just a single non-zero component of either **A** or **F**. In cases of more complicated geometry the problem can be reduced to solving multiple simple problems by means of mode decomposition.

In equation A1 we show the characteristic type of equation to be solved in order to calculate the non-zero component of the vector potentials, here in the case of a Z directed magnetic source.

$$F(x, y, z) \propto \hat{z} \iint_{-\infty}^{\infty} f(x, y, z, k_x^2 + k_y^2) e^{-i(k_x x + k_y y)} dk_x dk_y \quad (A1)$$

where $\hat{z} = i\mu\omega$, $\mu$ is the magnetic permeability, $\omega$ is the angular frequency, $k_x$ and $k_y$ is the wave number in the x and y direction.
We recognize this double integral as an inverse 2D Fourier transform, which can be reduced to a single integral by applying an inverse Hankel transform instead. This is possible since $f$ is cylinder symmetric given its dependence on $k_x^2 + k_y^2$ and allows for the equation to be rewritten in simpler form:

$$F(\rho, z) \propto \hat{z} \int_0^\infty f(\lambda, z) \cdot \lambda \cdot J_0(\lambda\rho) d\lambda \quad (A2)$$

where $\lambda$ is the space frequency, and $J_0$ is the 0th order Bessel function. This type of inverse Hankel transform integral can be evaluated very efficiently using digital filter methods as described by Johansen and Sørensen (1979) and Christensen (1990). Given the ability to efficiently solve inverse Hankel transform type integrals allows for fast calculation of any EM field component, by simple means of differentiation of the vector potential equations. In the calculation of the type of integral in equation A2 the most time consuming part becomes the evaluation of the integrand function $f$. This function includes a frequency dependent reflection coefficient, essentially accounting for the amplitudes of up- and down- going damped waves within the layered halfspace. This is calculated from a recurrence relation going up through the layering which becomes the effective bottlenecking term, making the calculation time scale linearly with the number of layers in the half space. In the case of frequency domain data these are simply the transmitted frequencies of the system, whereas the frequency space needs to be more systematically sampled for time domain problems. In this latter case we solve for a data type specific number of logarithmically spaced frequencies per decade and transform into the time domain using an inverse sine/cosine digital filter transform (Johansen and Sørensen (1979), Newman, Hohmann and Anderson 1986). This approach is
much more efficient than a standard inverse fast Fourier transform and for late
times we also find it more stable and accurate than the faster Gaver-Stehfest
inverse Laplace transform (Knight and Raiche 1982). To get to the final time
domain result from the discretely transformed frequency domain data we apply
an interpolating bicubic spline to obtain results for the particular points in time
of interest.

For the accuracy of the resulting responses we note that this is
determined mainly by the sampling density of the digital filters. Higher density
implies more computations, meaning that a performance/accuracy compromise
has to be made. We allow for the user to manually tune this trade off, but provide
default settings for an error of around 0.3% for frequencies up to 100 – 200 kHz
and for a time domain error of less than 1% for times as early as 4 µs. These
frequencies and times are somehow dependent on the overall conductivity of
the ground and thus high frequencies and early times are more inaccurate for
very high resistivities then for low resistivities.

Ever since the codebase was started more than 15 years ago, all
responses and common routines have been routinely validated against both
analytical expressions and other forward modeling algorithms.

## Modeling the system transfer function for TEM systems

An important aspect of modeling a complete TEM system is accounting for the
transmitter waveform. In the case of a step function this is done directly by
integration in the frequency domain (Bracewell 1986), while the general case of
arbitrary waveform is handled by convolution in the time domain using the
approach of Fitterman and Anderson (1987). These authors show that an
arbitrary waveform can be accounted for by a simple rewrite of the convolution
so it is expressed as a sum over a numerical differentiation of a piece wise
linear waveform.

Of equal importance is accounting for the receiver frequency characteristics.

Not only does a receiver coil itself have a frequency characteristic, but further
low-pass filters are typically also applied and both can have a significant
influence on the measured signal over the dynamic range of interest. Effersø,
Auken and Sørensen (1999) describes how low and high pass filtering are
implemented in the frequency domain before transforming to the time domain.

This is done by a simple multiplication of a filter response function, here in the
form of a Butterworth filter, with the same approach used for modeling of low
pass filters for time domain induced polarization (TDIP) responses. The effect
of finite gate width is included in the modeling by integrating cubic spline
functions over the width of the gate. If finite gates are not used the response at
the gate center times are calculated from a local interpolation using cubic
splines. Some instruments further utilize a special gate right after the receiver
coil. This gate serves to prevent strong primary fields from saturating the
receiver amplifiers during transmitter turn on and it is modeled by convolving a
shifted heavy side function with the step response.

Using SkyTEM as an example the full procedure for forward modeling
the response becomes: (1) the filter effect of the receiver coil is modeled in the
frequency domain as a simple product, (2) the frequency domain response is
transformed to the time domain using a Hankel transform, (3) the front gate is
convolved with the step response, (4) low pass filters in the receiver instrument
are applied as a new convolution, (5) a piece wise linear waveform is applied
607 by numerical differentiation and finally (6) gates are calculated by integrating

608 the response over the length of the gates.

609

610
APPENDIX B – LEAST SQUARES INVERSION

The inversion is performed iteratively, by following the established practice of linearized approximation of the non-linear forward mapping of the model to the data space, by the first term of the Taylor expansion (Menke 1989b, Auken et al. 2005a). The n+1th update of the model vector $\mathbf{m}_{n+1}$ is obtained by:

$$ \mathbf{m}_{n+1} = \mathbf{m}_n + [\mathbf{G}'_n \mathbf{C}'^{-1} \mathbf{G}'_n + \lambda_n \mathbf{I}]^{-1} \cdot [\mathbf{G}'_n \mathbf{C}'^{-1} \delta \mathbf{d}'_n] \quad (B1) $$

Where the Jacobian $\mathbf{G}'_n$, the data vector update $\delta \mathbf{d}'_n$ and the covariance matrix $\mathbf{C}'$ incorporate both the a priori and the roughness constraints and are defined as:

$$ \mathbf{G}'_n = \begin{bmatrix} \mathbf{G}_n \\ \mathbf{P} \\ \mathbf{R} \end{bmatrix} \quad (B2) $$

$$ \delta \mathbf{d}'_n = \begin{bmatrix} \delta \mathbf{d}_n \\ \delta \mathbf{m}_n \\ \delta \mathbf{r}_n \end{bmatrix} = \begin{bmatrix} \mathbf{d}_n - \mathbf{d}_{\text{obs}} \\ \mathbf{m}_n - \mathbf{m}_{\text{prior}} \\ -\mathbf{R} \mathbf{m}_n \end{bmatrix} \quad (B3) $$

$$ \mathbf{C}' = \begin{bmatrix} \mathbf{C}_{\text{obs}} & 0 & 0 \\ 0 & \mathbf{C}_{\text{prior}} & 0 \\ 0 & 0 & \mathbf{C}_{\mathbf{R}} \end{bmatrix} \quad (B4) $$

In equation B2 $\mathbf{G}_n$ represents the Jacobian of the forward mapping. For 1D LCI/SCI solutions the matrix $\mathbf{G}_n$ is computed by differentiation of forward responses $\mathbf{F}_{1D}: \mathbf{G}_n^{ij} = (\mathbf{F}_{1D}^i(\mathbf{m} + \Delta\mathbf{m}_j) - \mathbf{F}_{1D}^i(\mathbf{m}))/(\Delta\mathbf{m}_j)$, for datum $i$ and model parameter $j$. On the contrary, in the 2D DC/IP implementation the Jacoabin is computed with the adjoint method approach (Fiandaca et al. 2013b). $\mathbf{P}$ is the matrix with dimension $N_m \times N_m$ ($N_m$ being the number of model parameters), necessary to impose the constraints on the a priori values. $\mathbf{R}$ is the roughness matrix, in which each row represents one roughness constraint
(each row is zero everywhere except for the two elements corresponding to constrained parameters, elements being equal to 1 and -1).

The constraints also appear as extra terms in the definition of the data vector update $\delta d'_n$ (equation B3), that is composed of the distance $\delta d_n$ between the n-th forward response $d_n$ and the observed data $d_{obs}$, the distance $\delta m_n$ between the n-th model vector $m_n$ and the a priori model vector $m_{prior}$ (used also as starting model for the iterative procedure) and the roughness of the n-th model vector $\delta r_n = -Rm_n$.

The covariance matrix $C'$ of equation B4 is defined in terms of the covariance on the observed data $C_{obs}$, the covariance on the a priori information $C_{prior}$ and the covariance on the roughness constraints $C_R$. All three matrices are considered diagonal; the elements of $C_{prior}$ and $C_R$ control the strength of the constraints, while the elements of $C_{obs}$ reflect the noise content of the data.

The matrices $P$, $R$, $C_{prior}$ and $C_R$, as well as the vectors $\delta m_n$ and $\delta r_n$ are split in two parts in order to account for prior/roughness constraints on both the primary parameters (e.g. resistivity and thickness) and on depths. In fact, the depths are not included directly in the model space, and a special formulation of the matrices is needed to include the depth prior/roughness constraints in equation B1 (see Auken et al. 2005a for details).

In equation B1 the parameter $\lambda_n$ is the Marquart damping parameter (Marquart 1963), iteratively updated to stabilize the inversion process through an adaptive algorithm that damps the step size. This algorithm uses as damping value $\lambda_n$ the maximum diagonal value of the $G_n^{T}C^{-1}G'_n$ matrix, reduced by an iterative-dependent scaling factor $f_n$.

$$ \lambda_n = \max \operatorname{diag}(G_n^{T}C^{-1}G'_n) \cdot f_n \quad (B5) $$
The damping is used to stabilize the solution of the linear system of equations (B1), but a constraint on the “step size” of the iteration, i.e. the size of the model update ||m<sub>n+1</sub> - m<sub>n</sub>||, is also imposed through f<sub>n</sub> (the step size is reduced by increasing f<sub>n</sub>). When increasing the iteration number n, the scaling factor f<sub>n</sub> is decreased and the step size is increased ensuring a robust and efficient damping scheme.

The object function minimized by equation B1 is expressed by

$$Q = \left( \frac{[\delta d^T C d^{-1} \delta d^T]}{N_d + N_m + N_R} \right)^{\frac{1}{2}} \qquad (B6)$$

In which N<sub>d</sub>, N<sub>m</sub>, N<sub>R</sub> represent the number of data points, the number of model parameters and the number of constraints: The output models are then balanced between the data (through the forward response, i.e. the physics), the a priori constraints and the roughness constraints.

Finally the covariance of the estimator error C<sub>est</sub> (Tarantola and Valette 1982b) is used in AarhusInv to estimate the resolution of the inverted model by using its expression for linear mappings on the last iteration N<sub>ite</sub>:

$$C_{est} = (G'_{N_{ite}} C^{-1} G'_{N_{ite}})^{-1} \qquad (B7)$$
677

678

679 Reference List

680

681 Auken E. and Christiansen A.V. 2004. Layered and laterally constrained 2D

682 inversion of resistivity data. *Geophysics* **69**, 752-761.

683 Auken E., Christiansen A.V., Jacobsen B.H., Foged N., and Sørensen K.I.

684 2005a. Piecewise 1D Laterally Constrained Inversion of resistivity data.

685 *Geophysical Prospecting* **53**, 497-506.

686 Auken E., Christiansen A.V., Jacobsen L., and Sørensen K.I. 2008. A

687 resolution study of buried valleys using laterally constrained inversion of TEM

688 data. *Journal of Applied Geophysics* **65**, 10-20.

689 Auken E., Christiansen A.V., Viezzoli A., Fitzpatrick A., Cahill K., Munday T.,

690 and Berens V. 2009a. Investigation On The Groundwater Resources Of The

691 South Eyre Peninsula, South Australia, Determined From Laterally

692 Constrained Inversion Of Tempest Data. ASEG 2009 - Adelaide.

693 Auken E., Christiansen A.V., Viezzoli A., Fitzpatrick A., and Munday T. 2009b.

694 Laterally constrained inversion of TEMPEST data from Eyre Peninsula area,

695 South Australia. EAGE Near surface 2009.

696 Auken E., Christiansen A.V., Westergaard J.A., Kirkegaard C., Foged N., and

697 Viezzoli A. 2009c. An integrated processing scheme for high-resolution

698 airborne electromagnetic surveys, the SkyTEM system. *Exploration*

699 *Geophysics* **40**, 184-192.

700 Auken E., Nebel L., Sørensen K.I., Breiner M., Pellerin L., and Christensen

701 N.B. 2002. EMMA - A Geophysical Training and Education Tool for

702 Electromagnetic Modeling and Analysis. *Journal of Environmental &

703 *Engineering Geophysics* **7**, 57-68.

704 Balch S.J., Boyko W.P., and Paterson N.R. 2003. The AeroTEM airborne

705 electromagnetic system. *The Leading Edge* **22**, 562-566.

706 Behroozmand A.A., Auken E., Fiandaca G., and Christiansen A.V. 2012a.

707 Improvement in MRS parameter estimation by joint and laterally constrained

708 inversion of MRS and TEM data. *Geophysics* **74**, WB191-WB200.

709 Behroozmand A.A., Auken E., Fiandaca G., Christiansen A.V., and

710 Christensen N.B. 2012b. Efficient full decay inversion of MRS data with a

711 stretched-exponential approximation of the T2* distribution. *Geophysical*

712 *Journal International* **190**, 900-912.

713 Bracewell R.N. 1986. *The Fourier transform and its applications*, 2

714 edn.McGraw-Hill, Inc.

715 Brodie R. and Sambridge M. 2006a. A holistic approach to inversion of

716 frequency-domain airborne EM data. *Geophysics* **71**, G301-G312.
717 Brodie R. and Sambridge M. 2006b. A holistic approach to inversion of time-domain airborne EM. ASEG Extended Abstracts, 1-4.

719 Christensen N.B. 1990. Optimized fast Hankel transform filters. Geophysical Prospecting 38, 545-568.

721 Christensen N.B. 2002. A generic 1-D imaging method for transient electromagnetic data. Geophysics 67, 438-447.

723 Christensen N.B. and Auken E. 1992. Simultaneous electromagnetic layered model analysis. Proceedings of the interdisciplinary inversion workshop 1, Geoskrifter, 49-56.

726 Christiansen A.V. and Auken E. 2012. A global measure for depth of investigation. Geophysics 77, WB171-WB177.

728 Christiansen A.V., Auken E., and Viezzoli A. 2011. Quantification of modeling errors in airborne TEM caused by inaccurate system description. Geophysics 76, F43-F52.

731 Cox L.H., Wilson G.A., and Zhdanov M.S. 2010. 3D inversion of airborne electromagnetic data using a moving footprint. Exploration Geophysics 41, 250-259.

734 d'Ozouville N., Auken E., Sørensen K.I., Violette S., Marsily G.d., Deffontaines B., and Merlen G. 2008. Extensive perched aquifer and structural implications revealed by 3D resistivity mapping in a Galapagos volcano. Earth and Planetary Science Letters 2008, 518-522.

738 Danielsen J.E., Auken E., Jørgensen F., Søndergaard V.H., and Sørensen K.I. 2003. The application of the transient electromagnetic method in hydrogeophysical surveys. Journal of Applied Geophysics 53, 181-198.

741 Effersø F., Auken E., and Sørensen K.I. 1999. Inversion of band-limited TEM responses. Geophysical Prospecting 47, 551-564.

743 Fiandaca G., Auken E., Gazoty A., and Christiansen A.V. 2012a. Time-domain induced polarization: Full-decay forward modeling and 1D laterally constrained inversion of Cole-Cole parameters. Geophysics 77, E213-E225.

746 Fiandaca G., Ramm J., Binley A., Gazoty A., Christiansen A.V., and Auken E. 2013b. Resolving spectral information from time domain induced polarization data through 2-D inversion. Geophysical Journal International 192, 631-646.

749 Fitterman D.V. and Anderson W.L. 1987. Effect of transmitter turn-off time on transient soundings. Geoexploration 24, 131-146.

751 Gazoty A., Fiandaca G., Pedersen J., Auken E., and Christiansen A.V. 2012. Mapping of landfills using time-domain spectral induced polarization data: the Eskelund case study. Near Surface Geophysics 10, 563-574.
754 Gribenko A., Wan L., and Zhdanov M.S. 2010. Three-dimensional inversion of magnetotelluric data for complex resistivity. Proceeding of the CEMI Annual Meeting, 269-294.

757 Günther T., Liebau J., Akca I., and Mueller-Petke M. 2011. Block joint inversion (BJI) of MRS and DC resistivity soundings for aquifer imaging at the North sea island Borkum. Leicester.

760 Günther T., Rücker C., and Spitzer K. 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography - II. Inversion. *Geophysical Journal International* **166**, 506-517.

763 Hertrich M., Braun M., and Yaramanci U. 2005. Magnetic resonance soundings with separated transmitter and receiver loops. *Near Surface Geophysics* **3**, 141-154.

766 Johansen H.K. and Sørensen K.I. 1979. Fast Hankel transforms. *Geophysical Prospecting* **27**, 876-901.

768 Jørgensen F. and Sandersen P.B.E. 2006. Buried and open tunnel valleys in Denmark-erosion beneath multiple ice sheets. *Quaternary Science Reviews* **25**, 1339-1363.

771 Jørgensen F., Scheer W., Thomsen S., Sonnenborg T.O., Hinsby K., Wiederhold H., Schamper C., Roth.B., Kirsch R., and Auken E. 2012. Transboundary geophysical mapping of geological elements and salinity distribution critical for the assessment of future sea water intrusion in response to sea level rise. *Hydrology and Earth System Sciences* **16**, 1845-1962.

776 Kemna A., Binley A., Ramirez A., and Daily W. 2000. Complex resistivity tomography for environmental applications. *Chemical Engineering Journal* **77**, 11-18.

779 Kirkegaard C. and Auken E. 2014. A parallel, scalable and memory efficient inversion code for very large scale airborne EM surveys. *Geophysical Prospecting Accepted*.

782 Knight J.H. and Raiche A.P. 1982. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method. *Geophysics* **47**, 47-50.

785 Lane R., Green A.G., Golding C., Owers M., Pik P., Plunkett C., Sattel D., and Thorn B. 2000. An example of 3D conductivity mapping using the TEMPEST airborne electromagnetic system. *Exploration Geophysics* **31**, 162-172.

788 Loke M.H. 2011. RES3DINV, www.geoelectrical.com.

789 Loke M.H. and Barker R.D. 1996. Rapid least squares inversion of apparent resistivity pseudosections by a quasi-Newton method. *Geophysical Prospecting* **44**, 131-152.
792 Loke M.H., Chambers J.E., and Ogilvy R.D. 2006. Inversion of 2D spectral
793 induced polarization imaging data. *Geophysical Prospecting* **54**, 287-301.

794 Macnae J.C., Smith R., Polzer B.D., Lamontagne Y., and Klinkert P.S. 1991.
795 Conductivity-depth imaging of airborne electromagnetic step-response data.
796 *Geophysics* **56**, 102-114.

797 Marquart D. 1963. An Algorithm for Least Squares Estimation of Nonlinear
798 Parameters. *SIAM, Journal of Applied Mathematics* **11**, 431-441.

799 Menke W. 1989a. *Geophysical data analysis discrete inverse theory*. Academic Press.

801 Menke W. 1989b. *Geophysical data analysis: discrete inverse theory*. Academic Press.

803 Mohnke O. and Yaramanci U. 2002. Smooth and block inversion of surface
804 NMR amplitudes and decay times using simulated annealing. *Journal of
805 Applied Geophysics* **50**, 163-177.

806 Møller I., Verner H., Søndergaard V.H., Flemming J., Auken E., and
807 Christiansen A.V. 2009. Integrated management and utilization of
808 hydrogeophysical data on a national scale. *Near Surface Geophysics* **7**, 647-
809 659.

810 Mueller-Petke M. and Yaramanci U. 2010. QT inversion - Comprehensive use
811 of the complete surface NMR data set. *Geophysics* **75**, WA199-WA209.

812 Newman G.A. and Alumbaugh D.L. 1997. Three-dimensional massively
813 parallel electromagnetic inversion-I. Theory. *Geophysical Journal International*
814 **128**, 345-354.

815 Newman G.A., Hohmann G.W., and Anderson W.L. 1986. Transient
816 electromagnetic response of a three-dimensional body in a layered earth.
817 *Geophysics* **51**, 1608-1627.

818 Oldenburg D.W., Haber E., and Shekhtman R. 2013. Three dimensional
819 inversion of multisource time domain electromagnetic data. *Geophysics* **78**,
820 E47-E57.

821 Oldenburg D.W. and Jones F.H.M. 2011a. DCIP3D,
822 www.eos.ubc.ca/ubcgif/iag/.

823 Oldenburg D.W. and Jones F.H.M. 2011b. EM1DFM,
824 www.eos.ubc.ca/ubcgif/iag/.

825 Oldenburg D.W. and Jones F.H.M. 2011c. EM1DTM,
826 www.eos.ubc.ca/ubcgif/iag/.

827 Oldenburg D.W. and Li Y. 1994. Inversion of induced polarization data.
828 *Geophysics* **59**, 1327-1341.
829 Pelton W.H., Ward S.H., Hallof P.G., Sill W.R., and Nelson P.H. 1978. Mineral discrimination and removal of inductive coupling with multifrequency induced-polarization. *Geophysics* **43**, 588-609.

832 Podgorski J.E., Auken E., Schamper C., Christiansen A.V., Kalscheuer T., and Green A.G. 2013. Processing and inversion of commercial helicopter time-domain electromagnetic data for environmental assessments and geologic and hydrologic mapping. *Geophysics* **78**, E149-E159.

836 Reid J.E., Munday T., and Fitzpatrick A. 2007. High-Resolution Airborne Electromagnetic Surveying for Dryland Salinity Management: The Toolibin Lake SkyTEM Case Study, W.A. 1-5.

839 Rücker C., Günther T., and Spitzer K. 2006. Three-dimensional modelling and inversion of dc resistivity data incorporating topography - I. Modelling. *Geophysical Journal International* **166**, 495-505.

842 Santos F.A.M. 2004. 1-D laterally constrained inversion of EM34 profiling data. *Journal of Applied Geophysics* **56**, 123-134.

844 Sato H.K. 2000. Potential field from a dc current source arbitrarily located in a nonuniform layered medium. *Geophysics* **65**, 1726-1732.

846 Schamper C., Auken E., and Kirkegaard C. 2013. 3D Sheets Inversion With Accurate Modeling of AEM Systems. SAGA AEM, South Africa, 1-2.

848 Schamper C., Auken E., and Sørensen K.I. 2013. Coil response inversion for very early time modeling of helicopter-borne time-domain EM data and mapping of near-surface geological layers. *Geophysical Prospecting*.

851 Schamper C., Auken E., and Sørensen K.I. 2014. Coil response inversion for very early time modeling of helicopter-borne time-domain electromagnetic data and mapping of near-surface geological layers. *Geophysical Prospecting* DOI: 10.1111/1365-2478.12104.

855 Sengpiel K.P. and Siemon B. 2000. Advanced inversion methods for airborne electromagnetic exploration. *Geophysics* **65**, 1983-1992.

857 Siemon B., Stuntebeck C., Sengpiel K.P., Röttger B., Rehli H., and Eberle D.G. 2002. Investigation of Hazardous Waste Sites and their Environment Using BGR Helicopter-Borne Geophysical System. *Journal of Environmental & Engineering Geophysics* **7**, 169-181.

861 Sørensen K.I. and Auken E. 2004. SkyTEM - A new high-resolution helicopter transient electromagnetic system. *Exploration Geophysics* **35**, 191-199.

863 Sørensen K.I., Auken E., Christensen N.B., and Pellerin L. 2005. An Integrated Approach for Hydrogeophysical Investigations. New Technologies and a Case History. In: *Near-surface Geophysics Part II* pp. 585-603. SEG.
866 Sørensen K.I. and Larsen F. 1999. Ellog Auger Drilling: 3-in-one Method for
867 Hydrogeological Data Collection. *Ground Water Monitoring & Remediation* **19**,
868 97-101.

869 Tarantola A. and Valette B. 1982a. Generalized nonlinear inverse problems
870 solved using a least squares criterion. *Reviews of Geophysics and Space
871 Physics* **20**, 219-232.

872 Tarantola A. and Valette B. 1982b. Inverse problems = quest for information.
873 *J.of Geophysics* **50**, 159-170.

874 Tartaras E., Zhdanov M.S., Wada K., Saito A., and Hara T. 2000. Fast
875 Imaging of TDEM data based on S-inversion. *Journal of Applied Geophysics*
876 **43**, 15-32.

877 Thomsen R., Søndergaard V.H., and Sørensen K.I. 2004. Hydrogeological
878 mapping as a basis for establishing site-specific groundwater protection zones
879 in Denmark. *Hydrogeology* **12**, 550-562.

880 Viezzoli A., Christiansen A.V., Auken E., and Sørensen K.I. 2008. Quasi-3D
881 modeling of airborne TEM data by Spatially Constrained Inversion.
882 *Geophysics* **73**, F105-F113.

883 Viezzoli A., Munday T., Auken E., and Christiansen A.V. 2010. Accurate quasi
884 3D versus practical full 3D inversion of AEM data - the Bookpurnong case
885 study. *Preview* **149**, 23-31.

886 Vignoli G., Fiandaca G., Christiansen A.V., Kirkegaard C., and Auken E. 2013.
887 Sharp Spatially Constrained Inversion (sSCI) with applications to transient
888 electromagnetic data. *Submitted to Geophysical Prospecting*.

889 Wannamaker P.E., Hohmann G.W., and SanFilipo W.A. 1984.
890 Electromagnetic modeling of three-dimensional bodies in layered earths using
891 integral equations. *Geophysics* **49**, 60-74.

892 Ward S.H. and Hohmann G.W. 1988. Electromagnetic theory for geophysical
893 applications. In: *Electromagnetic methods in applied geophysics*, Vol. 1 (ed.
894 M.N. Nabighian), pp. 131-311. Society of exploration geophysicists (SEG).

895 Weichman P.B., Lovely E.M., and Ritzwoller M.H. 2000. Theory of surface
896 nuclear magnetic resonance with applications to geophysical imaging
897 problems. *Physical Review E - Statistical Physics, Plasmas, Fluids, and
898 Related Interdisciplinary Topics* **62**, 1290-1312.

899 Weidelt P. 1983. The harmonic and transient electromagnetic response of a
900 thin dipping dike. *Geophysics* **48**, 934-952.

901 Wisén R. and Christiansen A.V. 2005. Laterally and Mutually Constrained
902 Inversion of Surface Wave Seismic Data and Resistivity Data. *Journal of
903 Environmental & Engineering Geophysics* **10**, 251-262.
904 Witherly K.E., Irvine R.J., and Morrison E. 2004. The Geotech VTEM time
domain helicopter EM system. *ASEG Extended Abstracts* 1-4.

906 Xiong Z. 1989. Electromagnetic fields of electric dipoles embedded in a
907 stratified anisotropic earth. *Geophysics* **54**, 1643-1646.

908 Yoshioka K. and Zhdanov M.S. 2005. Three-dimensional nonlinear
909 regularized inversion of the induced polarization data based on the Cole-Cole
910 model. *Physics of the Earth and Planetary Interiors* **150**, 29-43.

911 Zhdanov M.S., Fang S., and Hursán G. 2000. Electromagnetic inversion using
912 quasi-linear approximation. *Geophysics* **65**, 1501-1513.

913 Zhdanov M.S. and Keller G.V. 1994. *The geoelectrical methods in
914 geophysical exploration*, 1 edn. Elsevier.

915 Zhou Q. 1989. Audio-frequency Electromagnetic Tomography for Reservoir
916 Evaluation. PhD dissertation, Lawrence Berkeley Laboratory, University of
917 California.

918
919
920
921
922
923 **Tables**

924
<table><thead><tr><th>Source group</th><th>Source type</th><th>Receiver type</th><th>Source position</th><th>Receiver position</th><th>Domain</th></tr></thead><tbody><tr><td rowspan="4">Dipole<br>Reference 1,2,3</td><td>Vertical magnetic dipole</td><td>Dipole: Any E and H field component</td><td>Anywhere</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Vertical electric dipole</td><td>Dipole: Any E and H field component</td><td>Anywhere</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Horizontal electric dipole</td><td>Dipole: Any E and H field component</td><td>Ground surface or in the air</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Horizontal magnetic dipole</td><td>Dipole: Any E and H field component</td><td>Ground surface or in the air</td><td>Ground surface or in the air</td><td>Frequency and time</td></tr><tr><td rowspan="5">Loop<br>Reference 1,2,3</td><td>Rectangular loop</td><td>Dipole: Any E and H field component</td><td>Ground surface or in the air</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Arbitrary segmented loop</td><td>Dipole: Any E and H field component</td><td>Ground surface or in the air</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Circular center</td><td>Dipole: Hz</td><td>Ground surface or in the air</td><td>Anywhere on the axis</td><td>Frequency and time</td></tr><tr><td>Circular offset</td><td>Dipole: Ex, Ey, Hx, Hy, Hz</td><td>Ground surface or in the air</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Circular in-loop</td><td>Hz</td><td>Ground surface or in the air</td><td>Ground surface or in the air</td><td>Frequency and time</td></tr><tr><td rowspan="4">Wire<br>Reference 1,2,3</td><td>Finite length grounded x-directed dipole</td><td>Dipole: Ex, Ez, Hy, Hz (or B)</td><td>Ground surface</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Finite length grounded y-directed dipole</td><td>Dipole: Ey, Ez, Hx, Hz (or B)</td><td>Ground surface</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Infinite x-directed line source</td><td>Dipole: Ex, Hy, Hz (or B)</td><td>Ground surface</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td>Infinite x-directed line source</td><td>Dipole: Ey, Hx, Hz (or B)</td><td>Ground surface</td><td>Anywhere</td><td>Frequency and time</td></tr><tr><td rowspan="3">DC<br>Reference 4,5, 6</td><td>Schlumberger</td><td>Apparent resistivity or potentials</td><td>Ground surface</td><td>Ground surface</td><td></td></tr><tr><td>Wenner</td><td>Apparent resistivity or potentials</td><td>Ground surface</td><td>Ground surface</td><td></td></tr><tr><td>Pole-pole</td><td>Apparent resistivity or potentials</td><td>Ground surface</td><td>Ground surface</td><td></td></tr></tbody></table>
<table><tbody><tr><td rowspan="3"></td><td>Arbitrary quad pole</td><td>Apparent resistivity or potentials</td><td>Ground surface</td><td>Ground surface</td><td></td></tr><tr><td>Arbitrary quad pole in the ground</td><td>Apparent resistivity or potentials</td><td>Electrodes at the surface or in the ground</td><td>Electrodes at the surface or in the ground</td><td></td></tr><tr><td>Arbitrary quad pole in cylinder symmetry</td><td>Apparent resistivity or potentials</td><td>Electrodes in the ground</td><td>Electrodes in the ground</td><td></td></tr><tr><td rowspan="2">DC/IP<br>Reference 7, 8</td><td>Arbitrary quad pole</td><td>Apparent resistivity/char geability or potentials</td><td>Ground surface</td><td>Ground surface</td><td>Time domain</td></tr><tr><td>Arbitrary quad pole in the ground</td><td>Apparent resistivity/char geability or potentials</td><td>Electrodes at the surface or in the ground</td><td>Electrodes at the surface or in the ground</td><td>Time domain</td></tr><tr><td>Magnetol-luric<br>Reference 1</td><td></td><td>Apparent resistivity or dipole Ex, Ey, Hx, Hy</td><td></td><td>Ground surface or in the ground</td><td></td></tr><tr><td>MRS<br>Reference 12, 13</td><td>Arbitrary segmented loop</td><td></td><td>Ground surface</td><td>Ground surface</td><td></td></tr><tr><td>Surface wave dispersion<br>Reference 14</td><td></td><td></td><td>Ground surface</td><td>Ground surface</td><td></td></tr></tbody></table>

Supported 2D and 3D responses

<table><thead><tr><td>Source group</td><td>Source type</td><td>Receiver type</td><td>Source position</td><td>Receiver position</td><td>Domain</td></tr></thead><tbody><tr><td>DC/IP (2D)<br>Reference 7, 8</td><td>Arbitrary quad poles</td><td>Apparent resistivity/char geability or potentials</td><td>Ground surface</td><td>Ground surface</td><td>Time</td></tr><tr><td>Loop – thin sheets (3D)<br>Reference 9, 10, 11</td><td>Arbitrary segmented loop</td><td>Dipole: Any E and H field component</td><td>Ground surface or in the air</td><td>Anywhere</td><td>Frequency and time</td></tr></tbody></table>

Table 1: Overview of supported data types and configurations. All time domain responses can be calculated as step or impulse response or with an arbitrary waveform. Abbreviations: 1D, 2D and 3D are the dimensionality of the response, E is electric field, H is magnetic field, Ex, Ey, Ez, Hx, Hy, Hz where x,y and z are the principal components, DC is direct current (geoelectric), DC/IP is DC with time domain induced polarization, MRS is Magnetic Resonance Sounding. Key references for the basic electromagnetic calculations include

¹Ward and Hohmann (1988), ²Wannamaker, Hohmann and SanFilipo (1984), ³Xiong (1989).
932 DC/IP responses in 1D are described by ⁴Zhdanov and Keller (1994) and ⁵Fiandaca et al. (2012) with the special case of electrodes in the layered halfspace which follows ⁶Sato (2000).
933 2D DC/IP uses a generalization of the approach described by ⁷Kemna et al. (2000) and is reported by ⁸Fiandaca et al. (2013). Sheets follows the theory of ⁹Weidelt (1983) and the implementation by ¹⁰Zhou (1989) with an additional time transform as discussed in ¹¹ Schamper, Auken and Kirkegaard (2013). Magnetic resonance soundings (MRS) follows the equations in ¹²Weichman, Lavely and Ritzwoller (2000) and is reported in ¹³Behroozmand et al. (2012b). Surface wave seismics is reported in ¹⁴Wisén and Christiansen (2005).

940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957

Optional Inversion parameters

<table><thead><tr><td>System type</td><td>Additional Parameters</td><td>Purpose</td></tr></thead><tbody><tr><td>All airborne TEM system</td><td>Receiver altitude/elevation</td><td>Correction for inaccurate determination of altitude</td></tr><tr><td>All airborne TEM systems</td><td>Data shift</td><td>Correction of data level</td></tr><tr><td>Fixed wing TEM systems</td><td>Receiver pitch roll and position</td><td>Correction for bird pitch, roll and position</td></tr><tr><td>MRS</td><td>Water content, T<sub>2</sub>* and stretch parameter (T<sub>2</sub>* distribution)</td><td>Additional parameters in model space</td></tr></tbody></table>
<table><tr><td>DC/IP</td><td>Chargeability, c and time constant (Cole-Cole parameters)</td><td>Additional parameters in model space</td></tr><tr><td>Sheets</td><td>Sheet position, size, strike and dip</td><td>Additional parameters in model space</td></tr></table>

958

959

<caption>Table 2: Additional inversion model parameters for relevant system types. The base model parameterization includes resistivity, thickness and depth</caption>

960

961
962

<table><thead><tr><td>Layer</td><td>Lithology</td><td>Resistivity (ohmm)</td><td>Thickness (m)</td><td>Water content (m³/m³)</td><td>T2* (s)</td><td>C</td></tr></thead><tbody><tr><td>1</td><td>Till, some water</td><td>40</td><td>30</td><td>0.1</td><td>0.1</td><td>1</td></tr><tr><td>2</td><td>Dry sand</td><td>300</td><td>20</td><td>0.02</td><td>0.2</td><td>1</td></tr><tr><td>3</td><td>Saturated sand</td><td>80</td><td>20</td><td>0.4</td><td>0.2</td><td>1</td></tr><tr><td>4</td><td>Clay</td><td>5</td><td></td><td>0.4</td><td>0.02</td><td>1</td></tr></tbody></table>

963

964 Table 3: The synthetic model used for mixed modeling of MRS and airborne TEM data

965

966

967 Figure Captions

968 Figure 1: Illustration of the modular algorithm design. All modules are implemented in a general manner to support any data type, except for the data file input and forward response modules which contains data type specific branches.

971

972

973 Figure 2: Joint surface and borehole DC/IP LCI inversion result. The model sections on the left side are faded below the depth on investigation and the right side of the figure shows the corresponding standard deviation sections. The position and depth of the logged borehole included in the inversion is marked by a black square and results from 2 boreholes is superimposed on the left hand figures.

978

979 Figure 3: SkyTEM SCI inversion results. (a) shows the flight lines of the dataset, (b) is the depth of investigation, (c) the modeled resistivity at a depth of 74-88 meters and (d) is the model parameter analysis of (c).

982

983 Figure 4: Results of MRS and Airborne TEM step-wise inversion. The green lines shows the true model and the black lines are the inverted model results. The dashed blue lines indicate
985 68% confidence intervals, as obtained from model parameter analysis, with a red line
986 superimposed to indicate uncertainty so large that it goes out of the scale of the figure. We
987 note that for the resistivity the sensitivity analysis is based on TEM data alone.

988
989
990 Figure 5: MRS and Airborne TEM full joint inversion results. The green lines show the true
991 model parameters and the black lines the inverted parameters. The dashed blue lines indicate
992 68% confidence intervals, as obtained from model parameter analysis.

993
994 Figure 6: Thin sheet inversion results where the sheet is buried under a conductive
995 overburden. (a) 3D view of the thin sheets with the starting, true and final estimated models
996 (those 2 last ones are confounded), the red dots corresponds to the airborne TEM sounding
997 locations; (b) 2D top view of the same sheet models; (c) Data profile with the synthetic
998 generated data (noise has been added at late times), the response of the starting model, and
999 the response of the final model.

1000
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>1016 Figure 6</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>