Three-Dimensional Inversion of Induced Polarization Effects in Airborne Time Domain Electromagnetic Data Using the GEMTIP Model

Abstract

This paper discusses the physical and mathematical principles of the airborne induced polarization (IP) method. The possibility of extracting information about the IP properties of rocks from airborne survey data has become a subject of active research recently. We introduce a method for the joint inversion of the airborne EM data into the electrical conductivity and IP parameters based on the generalized effective-medium theory of induced polarization (GEMTIP). We also present the results of the inversion of the airborne EM data collected over the Echum Project Area, in Northwestern Ontario, Canada, into 3D conductivity and chargeability models. Obtaining IP physical property models from an airborne geophysical survey may result in a paradigm change in mineral exploration by pulling more information and value from airborne EM surveys.

1. Introduction

Modern airborne surveys collect electromagnetic (EM) and total magnetic intensity (TMI) data simultaneously. Three-dimensional inversion of all collected data provides important information about the geology in the survey area and maximizes the value of the airborne survey. Traditionally, AEM data have been inverted for electrical conductivity only. This paper considers a more comprehensive approach, which allows us to rigorously invert the AEM data in full 3D for the electrical conductivity and induced polarization (IP) parameters of the subsurface. This both adds additional parameters for interpretation and improves the accuracy of the recovered conductivity.

It was observed long ago that the IP effect could be detected by using not only galvanic but also inductive sources and receivers [1,2]. The IP effect has also been manifested by the visible negative responses in the time domain EM data recorded with the EM induction system [3,4]. These observations have stimulated the growing interest in studying the IP phenomenon in airborne EM data. Smith and Klein [4] were the first to demonstrate that recovering the IP parameters from airborne EM data was possible. They also observed a correlation between IP inversion results produced from the data collected using the airborne platform with an induction system and the data measured using conventional ground-based galvanic sources and receivers. Since then, the concept of airborne IP has been discussed extensively in several publications [5,6,7,8]. Viezzoli et al. [9] and Kaminski and Viezzoli [10] conducted numerical modeling of the IP effect in airborne data. They also developed an inversion method based on the Cole–Cole model and traditional 1D layered-earth approximation of the geoelectrical section.

These publications demonstrated the potential of using airborne data to study the subsurface’s IP parameters. However, since 1D modeling and inversion do not adequately represent the actual 3D structure of geological formations, the electrical conductivity and IP parameters determined by 1D models may represent a distorted image of the subsurface, the extent of which will not be known until 3D modeling is undertaken [11].

Yoshioka and Zhdanov [12] and Xu and Zhdanov [13] developed full rigorous three-dimensional methods for the inversion of the induced polarization data based on the Cole–Cole model. Zhdanov [14] and Alfouzan et al. [15] introduced the 3D inversion method for arbitrary GEMTIP model parameters, with the Cole–Cole model being a special simplified case.

These recent developments are based on understanding the fundamental fact that the electrical resistivity of rocks is not necessarily positive and frequency-independent, and can be a complex function varying with frequency. This complex function is often called the resistivity relaxation model. It depends on a combination of structural and physical parameters of the model and the fluids filling the porous space. The IP effect can be calculated simultaneously with the electromagnetic induction response by applying Maxwell’s equations to the geoelectrical model with the given resistivity relaxation function.

Many different models have been proposed to replicate the complex resistivity spectra seen in real rocks. These include, for example, Debye dispersion [16], but the most prevalent is the Cole–Cole model [17]. This model was introduced to geophysics by Pelton [18,19,20]. A review of this method was conducted by Seigel et al. [21,22]. Indeed, the Cole–Cole model is most widely applied for modeling IP data [12,23,24,25,26]. However, it does not directly connect the mineralization and rock composition to the induced polarization spectrum.

In 2008, Zhdanov [27,28] introduced the GEMTIP model based on the generalized effective-medium theory of the IP phenomenon, which does allow a direct connection of the rock matrix to the complex resistivity spectrum. In these papers, Zhdanov showed that the Cole–Cole model could be derived analytically from the generalized effective-medium theory for electromagnetics. Indeed, the most common form of the Cole–Cole model is a specific solution to the GEMTIP model when conductive inclusions are spherical, and there is only one type or phase of a mineral that is giving rise to the observed induced polarization. In the GEMTIP model, the petrophysical properties of the rock can be directly modeled, including mineralization, porosity and fluid content, anisotropy, and multiple mineral phases. This was confirmed by Burtman [29] and Zhdanov et al. [14,30], who used QEMScan mineralogical testing along with laboratory measurements of the resistivity relaxation functions of the rock samples.

In this paper, we introduce a method of three-dimensional joint inversion of airborne EM data into the electrical conductivity and IP parameters of the subsurface. This method uses the integral equation (IE) method in combination with the GEMTIP resistivity relaxation model for the simultaneous numerical modeling of both EM induction and IP effects in the predicted data. The inversion is based on the re-weighted regularized conjugate gradient method [12,31], which ensures the robust convergence of the inverse iterative process. We also employ the concept of the moving sensitivity domain to perform inversion on the large-scale airborne survey data [32,33].

As an illustration of the practical effectiveness of the developed method, we present the results of the inversion of VTEM data collected over the Echum Project Area, 54 km ENE of Wawa in Northwestern Ontario, Canada [34]. We inverted the VTEM dB/dt data into 3D conductivity and chargeability models. Obtaining multiple 3D physical properties models from a single airborne geophysical survey maximizes the data value. Several drilling targets were developed based on the results of this 3D interpretation.

2. GEMTIP Resistivity Relaxation Model

We note in the Introduction that, in a general case, the effective resistivity of rocks is not necessarily a constant and real number, but rather is complex and may vary with frequency. This complex behavior of electric resistivity gives rise to the phenomenon known as the induced polarization (IP) effect. In this situation, there is no difference in what kind of sources are used for generating and measuring the EM field—galvanic or inductive. The complex resistivity will affect the data in either of these cases. This simple fact makes it clear that airborne EM data are also affected by the IP phenomenon, and the complex resistivity can be, in principle, extracted from the airborne data. However, one has to develop the proper processing and inversion algorithms to solve this problem practically. The first step in solving this problem requires selecting the proper resistivity relaxation model.

A general approach to constructing the resistivity relaxation model is based on rock physics and the description of the medium as a composite heterogeneous multiphase formation [14,27].

In the paper by Burtman et al. [29], we introduced, for simplicity, the frequency-dependent complex resistivity for a two-phase model with elliptical inclusions, described by the following formula:

$${\rho }_e={\rho }_0{\left\{1+\frac{f}{3}{\sum }_{\alpha =x,y,z}\frac{1}{{\gamma }_{\alpha }}\left[1-\frac{1}{1+s_{\alpha }{\left(i\omega \tau \right)}^{C_1}}\right]\right\}}^{-1},$$

(1)

where ${\rho }_0$ is the DC resistivity (Ohm-m); $\omega$ is the angular frequency (rad/sec); $\tau$ is the time parameter; and C is the relaxation parameter. The coefficients ${\gamma }_{\alpha }$ and $s_{\alpha }\left(\alpha =x,y,z\right)$ are the structural coefficients defined by geometrical characteristics of the ellipsoidal inclusions used to approximate the grains:

$$s_{\alpha }=\raisebox{1ex}{$r_{\alpha }$}\!\left/ \!\raisebox{-1ex}{$\overline{a}$}\right.$$

(2)

and $\overline{a}$ is an average value of the equatorial (a_x and a_y) and polar (a_z) radii of the ellipsoidal grains, i.e.,

$$\overline{a}=\frac{\left(a_x+a_y+a_z\right)}{3},$$

(3)

$$r_{\alpha }=2\frac{{\gamma }_{\alpha }}{{\lambda }_{\alpha }},$$

(4)

where ${\gamma }_{\alpha }$ and ${\lambda }_{\alpha }$ are the diagonal components of the volume and surface depolarization tensors described in Burtman et al. [29].

In the case of spherical inclusions of radius a, ${\gamma }_{\alpha }=1/3$, $r_{\alpha }=a$, and $s_{\alpha }=1$ for $\alpha =x,y,z$. Therefore, Formula (1) can be simplified as follows:

$${\rho }_e={\rho }_0{\left\{1+3f\left[1-\frac{1}{1+{\left(i\omega \tau \right)}^C}\right]\right\}}^{-1}.$$

(5)

After some algebra, we can transform Formula (5) into a form similar to the conventional Cole–Cole formula for the effective resistivity:

$$\rho \left(\omega \right)=\rho \left(1-\eta \left(1-\frac{1}{1+{\left(i\omega \tau \right)}^C}\right)\right),$$

(6)

where $\rho$ is the DC resistivity (Ohm-m); $\omega$ is the angular frequency (rad/sec); $\tau$ is the time parameter; $\eta$ is the intrinsic chargeability [35]; and C is the relaxation parameter. The dimensionless intrinsic chargeability, $\eta$, characterizes the intensity of the IP effect. It can be demonstrated that chargeability $\eta$ is approximately related to parameter $f$ according to the following formula [14]:

$$\eta \approx \frac{3f}{1+3f}.$$

(7)

Note that the inversion will be run with respect to the conductivity; therefore, it is more convenient to express the conductivity using Formula (5):

$$\sigma \left(\omega \right)=\sigma \left(1+3f\left(1-\frac{1}{1+{\left(i\omega \tau \right)}^C}\right)\right),$$

(8)

In this case, the anomalous conductivity, $\Delta \sigma$, is equal to

$$\Delta \sigma =\sigma \left(\omega \right)-{\sigma }_b,$$

(9)

where ${\sigma }_b$ is the background conductivity.

Thus

$$\Delta \sigma =\sigma f\left(\eta ,\tau ,C\right)-{\sigma }_b=\sigma \left(1+3f\left(1-\frac{1}{1+{\left(i\omega \tau \right)}^C}\right)\right)-{\sigma }_b,$$

(10)

where the function $f\left(\eta ,\tau ,C\right)$ is equal

$$f\left(\eta ,\tau ,C\right)=\left(1+3f\left(1-\frac{1}{1+{\left(i\omega \tau \right)}^C}\right)\right).$$

(11)

We use Expressions (10) and (11) in our algorithms for forward modeling and inversion of AEM data.

3. Modeling and Inversion of Airborne Time Domain EM Data for Complex Resistivity

With the goal of finding the most geologically reasonable model, which also replicates the measured electromagnetic data, we developed a comprehensive workflow to help reduce non-uniqueness, provide better stability, and avoid local minimums in 3D inversion. The first step is one-dimensional (1D) approximate inversion, which is used as quality control for the data and is also used to develop an approximate background conductivity model. For inversions including chargeability, the conductivity, chargeability, time constant, and relaxation coefficient are inverted simultaneously. One-dimensional inversion is faster than 3D inversion and suffers from less non-uniqueness, although this often means that the correct model cannot be obtained. The approximation for the inversion is that for each transmitter–receiver pair, the earth is laterally invariant to infinity for modeling and sensitivity calculations. While this accelerates computations, this approximation is not accurate when lateral variations are present in the earth—often the very areas that we wish to invert accurately. Additionally, for small offset systems such as VTEM, SkyTEM, and HeliTEM, the inline (X) component cannot be modeled and will be ignored when using 1D modeling and inversion.

The second step is full inversion, which involves considering all the underground geometry and can make use of both the X and Z components of the data. This results in more accurate 3D models, especially in areas with complex geometries and geology. However, this process requires more complex algorithms to be developed and applied. Advanced 3D modeling and inversion methods for interpreting AEM data can be found in various publications such as [32,33,36].

Furthermore, the inversion algorithms that have been developed consider both the electromagnetic induction and induced polarization effects during the unified inversion process. The principles of 3D induced polarization modeling and inversion are based on the generalized effective-medium theory of the induced polarization effect (GEMTIP), which was described earlier. Consequently, the conductivity, chargeability, and time constant recovered during inversion using the GEMTIP model provide information obtained from the airborne EM survey.

3.1. Modeling Airborne Time Domain EM Data by Integral Equation Method

Several numerical modeling techniques can be used for 3D computer simulation of the airborne time-domain electromagnetic (TDEM) data [14]. The most widely used among them are the finite difference (FD), finite element (FE), and integral equation (IE) methods. This paper applies the IE method for airborne TDEM data modeling and inversion because it has several practical advantages over the more popular FD or FE techniques. First, IE forward modeling requires the calculation of the Green’s tensors for the background conductivity model, which can be precomputed only once and saved for multiple uses on every iteration of an inversion, which speeds up the computation of the predicted data [31]. Second, IE forward modeling and inversion require the discretization of the domain of inversion only, while in the framework of the FD or FE methods, one must discretize the entire modeling domain (including the areas in the air as well). Third, in the framework of the moving sensitivity domain approach [33], the modeling domain of the IE method is selected within the sensitivity domain, which speeds up the computations and results in a relatively fast but rigorous inversion method.

We summarize below the basic principles of the IE method as applied to the airborne EM data. The EM field recorded at the receivers can be represented as a sum of the background EM field, E^b, H^b, and the anomalous EM field, E^a, H^a:

$$E=E^b+E^a,$$

(12)

$$H=H^b+H^a.$$

(13)

The anomalous electromagnetic field in the frequency domain is related to the electric current, $\mathit{j}\left(\omega \right)=\Delta \sigma \left(\omega \right)E\left(\omega \right)$, induced in the modeling domain, D, according to the following integral formulas [31]:

$$E^a\left(r_j;\omega \right)={\int \int \int }_D{\widehat{G}}_E\left(r_j|r;\omega \right)·\Delta \sigma \left(r_j;\omega \right)E\left(r_j;\omega \right)dv=G_E\left[\Delta \sigma E\right],$$

(14)

$$H^a\left(r_j;\omega \right)={\int \int \int }_D{\widehat{G}}_H\left(r_j|r;\omega \right)·\Delta \sigma \left(r_j;\omega \right)E\left(r_j;\omega \right)dv=G_H\left[\Delta \sigma E\right],$$

(15)

where r_j is the receiver position in some Cartesian coordinate system; ${\widehat{G}}_E\left(r_j|r;\omega \right)$ and ${\widehat{G}}_H\left(r_j|r;\omega \right)$ are the electric and magnetic Green’s tensors defined for an unbounded conductive medium with the background (horizontally layered) conductivity ${\sigma }_b$; $G_E$ and $G_H$ are the corresponding Green’s linear operators; and the domain D represents a volume with an arbitrarily varying conductivity, $\sigma \left(r_j;\omega \right)={\sigma }_b\left(r_j;\omega \right)+\Delta \sigma \left(r_j;\omega \right)$, within a domain D (where ${\sigma }_b\left(r_j;\omega \right)$ is a horizontally layered earth background conductivity).

One can write Equation (14) for the virtual receivers at point r′ located within the modeling domain, $r^{\prime }\in D$, as

$$E^a\left(r^{\prime };\omega \right)={\int \int \int }_D{\widehat{G}}_E\left(r^{\prime }|r;\omega \right)·\Delta \sigma \left(\mathrm{r};\omega \right)E\left(r_j;\omega \right)dv.$$

(16)

Substituting Expression (16) into Formula (12), we arrive at the integral equation with respect to the total electric field, $E\left(r^{\prime }\right)$,

$$E\left(r^{\prime };\omega \right)=E^b\left(r^{\prime };\omega \right)+{\int \int \int }_D{\widehat{G}}_E\left(r^{\prime }|r;\omega \right)·\Delta \sigma \left(\mathrm{r};\omega \right)E\left(r_j;\omega \right)dv.$$

(17)

Solving (17), one can find the total electric field inside the modeling domain D. Finally, we substitute this total electric field, the solution of Equation (17), into Equations (15) and (13) to calculate the magnetic field in the receivers. It was shown by [37] that one can transform (17) into a contraction operator equation, which can be effectively solved by different numerical methods. This paper uses this contraction integral equation formulation of Hursán and Zhdanov [37] in both forward modeling and inversion of the airborne data.

We should note that the anomalous conductivity in integral Equation (14) through (17) is defined by Formula (10); thus, this IE-based modeling approach automatically takes into account the IP effect.

For calculating the magnetic field in the time domain, we apply a cosine transform from the frequency to the time domain using the digital filter approach of [38]. In the case of the time-domain step response, $H_{step}\left(r,t\right)$, it is calculated as follows:

$$H_{step}\left(r,t\right)=-\frac{2{\mu }_0}{\pi }{\int }_0^{\infty }\frac{\mathrm{Im}\left(H_{step}\left(\omega ,t\right)\right)}{\omega }\cos \left(\omega t\right)d\omega ,$$

(18)

which is evaluated using a digital filter. In practice, the current waveform in the transmitter can be described by a known function, I(t). The magnetic fields in the receivers, $H\left(r,t\right)$, are then found as a convolution of the spline of the derivative of the transmitter current $\partial I\left(t\right)/\partial t$ with the step response:

$$H\left(r,t\right)=\frac{\partial I\left(t\right)}{\partial t}\ast H_{step}\left(r,t\right).$$

(19)

3.2. The Inverse Problem of the Airborne IP Method

The inverse problem of the airborne IP method can now be written in the classic form of the operator equation:

$$d=A\left(m\right),$$

(20)

where m is a vector formed by the anomalous conductivities, $\Delta \sigma \left(r,\omega \right)$, within the targeted domain, and d is a vector formed by the observed magnetic field data in the receivers, $H\left(r_j,t\right)$.

The inversion is based on minimization of the Tikhonov parametric functional, $P^{\alpha }\left(m\right)$, with the corresponding stabilizing functional, S(m) [39]:

$$P^{\alpha }\left(m\right)={\Vert W_d\left(A\left(m\right)-d\right)\Vert }_{L_2}^2+\alpha S\left(m\right),$$

(21)

where $W_d$ is the data weighting matrix and $\alpha$ is a regularization parameter (used to balance the misfit and stabilizer terms in Equation (21).

There are several possible choices for the stabilizer, S(m), including minimum norm, S_MN, minimum support, S_MS, and minimum gradient support, S_MGS, stabilizing functionals. These stabilizers are determined as follows [31]:

$$S_{MN}\left(m\right)={\Vert W_m\left(m-m_{apr}\right)\Vert }_{L_2}^2={\int }_V{\left[W_m\left(r;\omega \right)\left(\Delta \sigma \left(r;\omega \right)-\Delta {\sigma }_{apr}\left(r;\omega \right)\right)\right]}^{2}dv,$$

(22)

$$S_{MS}\left(m\right)={\int }_V\frac{{\left[\Delta \sigma \left(r;\omega \right)-\Delta {\sigma }_{apr}\left(r;\omega \right)\right]}^2}{{\left[\Delta \sigma \left(r;\omega \right)-\Delta {\sigma }_{apr}\left(r;\omega \right)\right]}^2+e^2}dv,$$

(23)

and

$$S_{MGS}\left(m\right)={\int }_V\frac{\nabla \Delta \sigma \left(r;\omega \right)·\nabla \Delta \sigma \left(r;\omega \right)}{\nabla \Delta \sigma \left(r;\omega \right)·\nabla \Delta \sigma \left(r;\omega \right)+e^2}dv,$$

(24)

where W_m is the diagonal weighting matrix of the model parameters with a diagonal formed by the values of the weighting function, and $w_m\left(r;\omega \right)$ is a focusing parameter determining the degree of sharpness of the inverse conductivity images.

The minimum norm stabilizer, S_MN, minimizes the variations of the anomalous conductivity, $\Delta \sigma \left(r;\omega \right)$, from the a priori model, $\Delta {\sigma }_{apr}\left(r;\omega \right)$. The minimum support stabilizer, S_MS, minimizes the volume occupied by the anomalous conductivity; while the minimum gradient support stabilizer, S_MGS, selects the inverse models with sharp boundaries between the formations with different electrical properties. Thus, by selecting the proper stabilizing functionals, the user may emphasize different properties of the inverse models. The minimum norm stabilizer usually results in relatively smooth distributions of conductivity, while the focusing stabilizers, S_MS and S_GMS, generate models with sharp boundaries. We use regularized inversion with focusing stabilization, as this recovers models with sharper boundaries and higher contrasts than the regularized inversion with smooth (e.g., minimum norm) stabilization.

By minimizing the parametric functional, $P^{\alpha }\left(m\right)$, we find the solution to the airborne IP inverse problem. A standard technique to find a minimum of $P^{\alpha }\left(m\right)$ is to apply a gradient-type minimization method [31]. We use the re-weighted regularized conjugate gradient method (RCGM) to find the solution of the inverse problem. The iterative process is terminated when the misfit reaches the given level of noise in the data, ${\epsilon }_0$:

$$\varphi \left(m_N\right)={\Vert r_N\Vert }^2\le {\epsilon }_0^2.$$

The appropriate selection of the data and model parameters weighting matrices is very important for the success of the inversion. The data weights for inversion are based on a two-part model, comprising absolute and relative errors. The absolute noise level considers the instrument noise floor and prevents small data values close to 0 from being overly important. The relative error level represents errors such as tilt and flight height errors. The estimated errors are listed for each survey individually. The inverse of these errors is used as data weights.

The computation of the model weighting matrix, W_m, is based on sensitivity analysis. In the current paper, we select W_m as the square root of the sensitivity matrix for the model in each iteration:

$$W_m^{\left(n\right)}=\sqrt{diag\left(F_{m_n}^{*}{F}_{m_n}\right)}.$$

As a result, we obtain a uniform sensitivity of the data to different model parameters [31].

We apply the adaptive regularization method. The regularization parameter $\alpha$ is updated in the process of the iterative inversion as follows:

$${\alpha }_n={\alpha }_1{q}^{n-1}; n=1,2,3\cdots ;0<q<1.$$

In order to avoid divergence, we begin an iteration from a value of ${\alpha }_1$, which can be obtained as a ratio of the misfit functional and the stabilizer for an initial model, then reduce ${\alpha }_n$ according to the last formula on each subsequent iteration and continuously iterate until the misfit condition is reached:

$$r_{n0}^w=\Vert r_{n0}^w\Vert =\Vert W_d\left[A\left(m_{{\alpha }_{n0}}\right)-d\right]\Vert /\Vert W_{d}d\Vert \le \delta ,$$

(25)

where $r_{n0}^w$ is the normalized weighted residual and $\delta$ is the relative level of noise in the weighted observed data.

3.3. Fréchet Derivative Calculation Using the Quasi-Born Approximation

The main step of the RRCG method is the calculation of the variation (or gradient) of the misfit functional. This step requires the computation of the Fréchet derivative of the observed data with respect to model parameters. We assume, as above, that the conductivity within a 3D geoelectrical model can be presented by the background (horizontally layered) conductivity ${\sigma }_b$, and an arbitrarily varying conductivity, $\sigma ={\sigma }_b+\Delta \sigma$, within a domain D.

In this model, the anomalous magnetic field is produced by the anomalous conductivity distribution $\Delta \sigma$, and it can be calculated according to Formula (15). Using this integral representation, we can express the corresponding Fréchet derivative, F_H, as follows:

$$F_H\left(r_j|r;\omega \right)={\left.\frac{\partial H\left(r_j\right)}{\partial \Delta \sigma \left(r\right)}\right|}_{\Delta \sigma }={\left.\frac{\partial H^{\Delta \sigma }\left(r_j\right)}{\partial \Delta \sigma \left(r\right)}\right|}_{\Delta \sigma }=\frac{\partial G_H\left[\Delta \sigma H\right]}{\partial \Delta \sigma }\approx {\widehat{G}}_H\left(r_j|r;\omega \right)·E\left(r;\omega \right).$$

(26)

We can treat the electric field E⁽ⁿ⁾(r), found for iteration number n, as the electric field in the above equations for a subsequent iteration (n + 1), E(r) = E⁽ⁿ⁾(r).

In this case, the Fréchet derivatives at iteration number n can be found by direct integration from the last equation involving the electric field E⁽ⁿ⁾(r) computed on the current iteration:

$$F_H\left(r_j|r;\omega \right)=\approx {\widehat{G}}_H\left(r_j|r;\omega \right)·E^{\left(n\right)}\left(r;\omega \right).$$

(27)

Note that the electric field $E\left(r;\omega \right)$ is calculated using the rigorous IE forward modeling method at each iteration step to compute the predicted data (EM field at the receivers). Therefore, no extra computation is required to find the electric field for the Fréchet derivative calculation.

In the case of time domain EM data, the Fréchet derivatives also need to be transformed into the time domain. This can be accomplished by means of the same technique as that for the magnetic field, i.e., via the Fourier transform using digital filtering techniques. For example, in the case of a causal step turnoff response, the Fréchet derivative in the time domain can be expressed as follows:

$$f_H\left(r_j|r;t\right)=-\frac{2}{\pi }{{\displaystyle \int }}_0^{\infty }\frac{\mathrm{Im}\left(F_H\left(r_j|r;\omega \right)\right)}{\omega }\cos \left(\omega t\right)d\omega .$$

(28)

As a result, our inversion technique, based on the IE method, requires just one forward modeling on every iteration step, while the conventional inversion scheme requires, as a rule, at least three forward modeling solutions per inversion iteration (one to compute the predicted data, another one to compute the gradient direction, and at least one for optimal calculation of the iteration step). This approach results in a very efficient inversion method.

3.4. Calculation of the Fréchet Derivatives with Respect to the GEMTIP Model Parameters

In the current paper, the resistivity relaxation model of the subsurface is represented by the GEMTIP model describing a two-phase composite medium. This GEMTIP model is characterized by four parameters (conductivity of the host medium, ${\sigma }_0$, fraction volume f, time constant, $\tau$, and relaxation coefficient, C). Therefore, the Fréchet derivatives with respect to each GEMTIP parameter are required to invert the IP data.

The Fréchet derivative of the magnetic field with respect to the conductivity of the host medium can be computed as follows:

$$\begin{gathered} F_H^{{\sigma }_0}\left(r_j|r;\omega \right)=F_H\left(r_j|r;\omega \right)·\frac{\partial {\sigma }_e\left(r;\omega \right)}{\partial {\sigma }_0}, \\ \frac{\partial {\sigma }_e\left(r;\omega \right)}{\partial {\sigma }_0}=1+\frac{f\left(r\right)}{3}\left\{1-\frac{1}{1+{\left[i\omega \tau \left(r\right)\right]}^{C\left(r\right)}}\right\} \end{gathered}$$

(29)

where $F_H\left(r_j|r;\omega \right)$ is the Fréchet derivative defined by Formula (27).

The Fréchet derivative of the magnetic field with respect to the fraction volume can be computed as follows:

$$\begin{gathered} F_H^f\left(r_j|r;\omega \right)=F_H\left(r_j|r;\omega \right)·\frac{\partial {\sigma }_e\left(r;\omega \right)}{\partial f}, \\ \frac{\partial {\sigma }_e\left(r;\omega \right)}{\partial f}=1+\frac{{\sigma }_0}{3}\left\{1-\frac{1}{1+{\left[i\omega \tau \left(r\right)\right]}^{C\left(r\right)}}\right\}, \end{gathered}$$

(30)

The Fréchet derivative of the magnetic field with respect to the time constant can be determined as follows:

$$\begin{gathered} F_H^{\tau }\left(r_j|r;\omega \right)=F_H\left(r_j|r;\omega \right)·\frac{\partial {\sigma }_e\left(r;\omega \right)}{\partial \tau }, \\ \frac{\partial {\sigma }_e\left(r;\omega \right)}{\partial \tau }=\frac{{\sigma }_{0}i\omega C\left(r\right){\left[i\omega \tau \left(r\right)\right]}^{C\left(r\right)-1}\left\{\frac{f\left(r\right)}{3}-\left(1+\frac{f\left(r\right)}{3}\right){\left[i\omega \tau \left(r\right)\right]}^{C\left(r\right)}\right\}}{{\left(1+{\left[i\omega \tau \left(r\right)\right]}^{C\left(r\right)}\right)}^2} \end{gathered}$$

(31)

The Fréchet derivative of the EM fields with respect to the relaxation coefficient is equal to:

$$\begin{gathered} F_H^C\left(r_j|r;\omega \right)=F_H\left(r_j|r;\omega \right)·\frac{\partial {\sigma }_e\left(r;\omega \right)}{\partial C}, \\ \frac{\partial {\sigma }_e\left(r;\omega \right)}{\partial C}=\frac{{\sigma }_{0}f\left(r\right){\left[i\omega \tau \left(r\right)\right]}^{C\left(r\right)}ln\left(i\omega \tau \left(r\right)\right)}{3{\left\{1+{\left[i\omega \tau \left(r\right)\right]}^{C\left(r\right)}\right\}}^2} \end{gathered}$$

(32)

In order to compute the Fréchet derivatives in the time domain, one can use the cosine transform given by Formula (28) above.

4. Case Study: Inversion of Airborne Data Acquired in Wawa, Ontario, Canada

4.1. Airborne Data Collection and Geological Background of the Survey Area

The area of investigation is over the Wawa Greenstone Belt. This area has been previously examined with several geophysical methods in a brief abstract [34]. The regional geology is 2.89 to 2.70-billion-year-old Precambrian rock which starts northeast of Lake Superior and extends to the western edge of the Kapuskasing Horst structural zone of migmatized rock. Mafic to ultramafic bodies of various ages intrude into the metavolcanic–metasedimentary belt. Common economic mineralization in this area includes gold, silver, zinc, copper, and iron mineralization. Diamondiferous kimberlites and lamproites are observed in the southeastern part of the Wawa Greenstone belt. An excellent description of the regional and detailed geology can be found in [40]. Gold properties are found around the northwest periphery of the same granite-granodiorite batholith that occurs along the east side of the property. Figure 1 below shows the location of the survey area.

The airborne survey was flown in February of 2021. The survey employed the Geotech Time Domain EM (VTEM^TM Plus) full receiver-waveform streamed data recorded system. Figure 2 shows the full VTEM survey flight path in red, plotted on a satellite image. The VTEM system had a 64 m terrain clearance. The system was concentric, with the X (inline) and Z (vertical) receivers being located near the center and in the same plane as the 26 m diameter horizontal transmitter loop. The receivers measured the time derivative of the magnetic field over 43 time gates from 21 microseconds to 8 ms. The inversion used data starting at 42 microseconds to 8 ms. The full waveform system was used, which measured the entire period of the instrument’s waveform and then was post-processed to an idealized decay. The TMI magnetic sensor was 10 m above the transmitter loop.

4.2. VTEM Inversion

First, a 1D inversion was performed (using the z component only) to determine a background and starting model for the 3D inversion. This background model included conductivity, volume fraction, time constant, and relaxation coefficient. The 3D inversion was performed on a grid with 50 × 50 m cells in the horizontal direction and 14 cells in the vertical direction, from 10 m-thick at the surface to 60 m-thick at a maximum depth of 500 m. The final 3D inversion was run on a node with 16 processors and 196 GB RAM. The Matlab parallel toolbox was used to leverage all processors. The final inversion took about 30 h. The full inversion domain was 196 cells × 86 cells × 14 cells for a total of 235,984 active cells in the inversion domain, or nearly 1 million free parameters when the 4 IP parameters were considered.

The inversion was run until the data were fit to the estimated noise level in the observed data. The errors were estimated by means of a two-part error model given by Equation (33):

$${\mathbf{\epsilon }}^T=\left|{\mathit{d}}^{obs}\right|\ast {\mathbf{\epsilon }}^P+{\mathbf{\epsilon }}^A$$

(33)

where ${\mathbf{\epsilon }}^T$ is a vector of the total error in each data point, $\left|{\mathit{d}}^{obs}\right|$ is the magnitude of each point in the observed data, ${\mathbf{\epsilon }}^P$ is the estimated percent error in the dataset, ${\mathbf{\epsilon }}^A$ and is the estimated absolute error level or error floor in the data set. The inversion was run until the normalized chi squared metric reached 1:

$${\chi }^2=\frac{1}{N_d}{\left({\mathit{r}}^{\mathit{T}}\mathit{r}\right)}^{\frac{1}{2}}$$

(34)

where r is the error normalized residual:

$$\mathit{r}=\frac{{\mathit{d}}^{obs}-{\mathit{d}}^{pred}}{{\mathbf{\epsilon }}^T}.$$

(35)

The percentage noise levels were set to 15% for both the x and z channels, and the absolute errors (noise floor) were set to 0.005 nT/Am⁴ for x and 0.001 nT/Am⁴ for z. For conductivity-only inversions, db_z/dt data points that were less than 1/3 the noise floor were rejected completely. If IP was included, then instead, db_z/dt data points that had a magnitude less than 1/3 the noise level were rejected, but large negative values were kept. Positive and negative db_x/dt above the noise level were always retained.

4.2.1. Conductivity Results

We will start by examining the inverse model of conductivity. In Figure 3, there is a horizontal slice of the 3D conductivity model recovered at a depth of 100 m below the surface. The survey area exhibits a conductivity range between approximately 1 S/m (1 Ohm-m) to 1 × 10⁻⁴ S/m (10,000 Ohm-m), with much of the terrain being highly resistive (>1000 Ohm-m). The image displays both conductive and resistant lineaments running northwest to southeast. Some mild line stripping is evident in the figures, particularly in the northwest. The most noticeable feature in the data is the MPD Zinc Copper mineralization area, which appears as a strong conductor in the western third of the image at 709,000 mE and 5,341,500 mE.

Figure 4 and Figure 5 show the MPD mineralized area in more detail. Figure 4 shows a data profile over this mineralization along with a conductivity section. Observed and predicted db_x/dt data are shown here. In this location, the x-component responds very well to the compact, steeply dipping conductive body. Figure 5 shows vertical planes of conductivity, plus an isosurface at a constant conductivity of 0.5 S/m. Figure 5 shows the conductive core of the body.

4.2.2. Chargeability Inversion Results

This section presents the airborne-derived 3D chargeability model. We convert the GEMTIP parameter volume fraction (f) into chargeability according to relation (7). The chargeability from the airborne survey results typically do not and should not necessarily match with ground-based IP measurements. We believe this is primarily because of the greatly different bandwidths of ground-based systems (~1/4 Hz) and airborne systems (30 Hz). This causes two-orders-of-magnitude differences in the time constants of the media that they image. However, the information presented in the two is complementary, and the AEM systems provide useful chargeability measurements.

Figure 5 shows a horizontal slice from the 3D airborne-derived chargeability at a depth of 150 m below the surface. The large conductors are shown as chargeability lows, largely because the coupled EM and IP responses lead the IP response to be covered in very conductive areas. As is common with AEM chargeability measurements in our experience, much of the chargeability corresponds to lake bottom sediments. The chargeability is likely restricted to the near surface under these locations, but is shown from the near surface to 150 m depth in the inversion results. This is due to a lack of resolution in the inversion for the induced polarization. However, the AEM chargeability is capable of imaging petrophysical properties and can be used for mapping.

Figure 6 and Figure 7 show chargeability overlain on satellite imagery and a geologic map focusing on the southeast part of the survey area. Figure 6 shows that most of the chargeability corresponds to lake bottom sediments, but areas such as around 711,500 mE and 5,341,000 mN (red outline) show a high chargeability that is not associated with a lake and is of interest. Note that this anomaly corresponds to a major fault and contact between the mafic volcanics and granodiorite (Figure 7). There is also an obvious low chargeability relative to the background depicted in dark blue to the southeast that warrants further investigation.

Chargeability at a depth of 150 m overlain on a geologic map is presented in Figure 8 with the legend shown in Figure 9. Figure 10 shows a short data profile over the chargeable anomaly circled in Figure 8, along with recovered chargeability and conductivity sections. The negative transient between 711,400 and 711,600 mE clearly demonstrates polarizable material. For a zero transmitter–receiver offset system such as VTEM, there is no conductivity structure that is frequency-independent, which causes negative data. Figure 11 shows the decay of a single sounding point at 711,500 mE in Figure 10. The decay shows a coherent negative transient, which fits well with the inversion. Figure 12 shows the same section as Figure 10, but with conductivity-only inversion. The data fit is very poor, and the conductivity section is much less coherent. One can see a dramatic difference in the inverse images produced by the inversions without IP (Figure 12) and with IP (Figure 10). This result illustrates the importance of considering the IP effect in the AEM data inversion.

A single observed and predicted transient from 711,500 mE from Figure 10 is presented in Figure 11. The blue circles show the observed data, the red pluses show the predicted data, and the black asterisks show data that are nulled due to low amplitude. Profile of data and conductivity from inversion without induced polarization included is presented in Figure 12. The top panel shows the observed data as solid lines and the predicted data as dotted lines. The black asterisks indicate data that are nulled because of their low amplitude. The data fit poorly. The bottom panel shows a vertical cross section of conductivity.

5. Conclusions

In this paper, we have considered the possibility of extracting information about the IP properties of the rocks from airborne EM survey data. The solution to this problem required the development of a method of joint 3D modeling and inversion of airborne EM data based on the GEMTIP resistivity relaxation model. Full 3D inversion of the airborne EM data considering both the EM induction and IP effects increases the effectiveness of airborne geophysical methods in mineral exploration both by adding additional information (chargeability) and improving the accuracy of the recovered conductivity.

We have presented the results of the application of this method to the airborne EM data collected over the Echum Project Area in Northwestern Ontario, Canada. We have demonstrated that the results of these inversions correlate well with the known geology in the area. Several examples of potential targets have been suggested based on our understanding of the area. We have identified chargeability anomalies from airborne EM data that warrant further investigation for gold and disseminated sulfide mineralization. Until drilling is complete, we will not know if the target anomalies are sulfides or another source. However, including IP in the inversion of airborne data greatly improved the data fit and provided more accurate conductivity sections. Additionally, the induced polarization anomalies are clearly geologic in origin, and if they are due to clays or other alterations, this still can be important information for exploration.

Finally, we conclude that it is essential to consider the IP effect in the AEM data inversion in order to produce the correct conductivity model, and the airborne IP model can provide important complementary information to ground-based surveys.