Detecting induced polarisation effects in time-domain data - a modelling study using stretched exponentials

...

1Introduction¶

Some earth materials are chargeable because they can store charge when an electric field is applied by an electromagnetic (EM) source. This is often called the induced polarization (IP) phenomenon. These materials can have different polarization mechanisms which results in different IP characteristics as a function of frequency. This can be translated into a complex conductivity model such as the Cole-Cole conductivity model Cole & Cole (1941)Pelton et al. (1978)Tarasov & Titov (2013):

where $\siginf$ is the conductivity at infinite frequency, $\eta_{cc}$ is the chargeability, $\tau_{cc}$ is the time constant (s), and $c_{cc}$ is the frequency exponent; the subscript CC indicates Cole-Cole.

IP surveys have been successfully conducted in a variety of geoscience applications. For mining, IP surveys are recognized as a principal geophysical technique for finding disseminated sulphides or porphyry deposits Fink et al. (1990). Non-metallic materials such as clays and ice can also generate IP signals Grimm & Stillman (2015)Leroy & Revil (2009); this makes IP a useful technique in a range of environmental applications. Grounded DC-IP surveys have been successfully used for both mining and environmental applications for the past decades. Airborne EM (AEM) systems can also detect IP signals. In particular time-domain AEM surveys using a coincident-loop system sometimes display a negative transient; this is a distinctive IP signature Weidelt (1982). Compared to EM signals, these negatives (IP signals) are much smaller in amplitude. Hence, for the initial AEM systems, it was not clear if the measured negatives were signals from chargeable materials or if they were simply noise generated by power lines or electric fences Smith & Klein (1996). With time however, instruments have improved and the validity of negative transients as signal has been firmly established Macnae (2016)Viezzoli et al. (2017). For instance, consistent negatives were recorded over the Tli Kwi Cho kimberlite deposit with three different AEM systems Kang et al. (2017). As the quality of instrumentation improves, it is expected that more IP signals will be measured in airborne data. This ability provides motivation for developing methodologies that can extract chargeability information from airborne IP data. Various approaches, including simple curve-fitting, 1D inversions, and 3D inversions have been developed and successfully applied to field examples Kratzer & Macnae (2012)Kwan et al. (2015)Hodges & Chen (2014)Kaminski & Viezzoli (2017)Kang & Oldenburg (2017). There is a significant enthusiasm for the potential use of the airborne IP techniques in a variety of applications (e.g. mining and groundwater). However, setting proper expectations about the technique, and recognizing its limitation based upon the physics and the current system specifications, is crucial because neither overselling nor underselling the technique is beneficial for the community.

Macnae (2016) investigated the physics of airborne IP and its practical aspects using a simple thin-sheet solution. A main conclusion from his study was that airborne IP is effectively a surficial clay mapper ($z<$ 100 m). Viezzoli et al. (2017), showed the potential that a deeper chargeable target, such as a mineral deposit ($z>$ 100 m), can be detected. That work however, was based upon analysis using 1D simulations. Hence, there is disagreement about the potential depth of investigation of the airborne IP technique. Although the approximate thin-sheet solution and semi-analytic 1D solutions used in Macnae (2016) and Viezzoli et al. (2017), respectively, can illustrate some meaningful concepts with respect to airborne IP, these approaches are limited in their ability to model the physics in the presence of complex conductivity structures. For instance, the finite size of the chargeable structure (e.g. width and length) is not taken into account in either approach. Investigating the feasibility of airborne IP in realistic geologic settings requires the use of 3D numerical simulations that solve the full Maxwell’s equations.

In this paper, we first develop a convolutional time-domain EM (TEM) simulation code using a stretched exponential (SE) conductivity function Kohlrausch (1854). We then use this code to investigate four main questions related to the feasibility of the airborne IP under ranges of circumstances:

• How does chargeable material in the subsurface generate negative transients in coincident loop systems?

• How does the background conductivity affect the IP signals?

• To what depth can we expect to detect a chargeable body?

• What are the characteristics of detectable chargeable materials in AEM data?

For our feasibility study, we limit our attention to detectability of IP signals, and we do not address issues of resolvability of chargeable structures in the inversion; that issue is beyond the scope of this study.

2Simulating airborne IP data using a streched exponential¶

With a complex conductivity, $\sigma(\omega)$, the current density, $\J$, in the frequency domain, can be written as:

where $\E$ is the electric field (V/m). In the time-domain, the current density, $\j$, is:

where $\otimes$ is a convolution. Then Maxwell’s equations can be written as

where $\b$ is the magnetic flux density (Wb/m²) and $\j_s$ (A/m²) is the current source; μ is the magnetic permeability (H/m). By discretizing and solving the above equations in 3D, we can compute TEM data that include IP effects Marchant et al. (2014)Marchant (2015). For the discretization of eqs. (2)-(5), an open-source geophysical simulation and inversion package, SimPEG, is used Cockett et al. (2015). The developed SimPEG-EMIP code works for both 3D tensor meshes and 2D/3D cylindrical meshes Heagy et al. (2017). Although not shown in the paper, the SimPEG-EMIP code can handle arbitrary waveforms such that user can input actual system waveforms used for their own case studies. For further details about solving the convolutional form of Maxwell’s equations, see Appendix A. The code is tested with an analytic solution described in Appendix B.

For a time-dependent conductivity, $\sigma(t)$, we use the stretched exponential (SE) model rather than the Cole-Cole model defined in the frequency-domain (eq. 1). The SE conductivity for a step-off function, $1-u_{step}(t)$, can be written as

where $u_{step}(t)$ is the Heaviside step function, $\sigma_0=\siginf(1-\eta_{se})$ is the DC conductivity, and subscript SE stands for stretched exponential. We want to obtain $\sigma_{se}$ from eq. 6. Taking the derivative with respect to time and multiplying by -1 yields:

where $\delta (t)$ is the Dirac-Delta function. Evaluating eq. (7) with eq. (6) results in

A main reason why we used the SE conductivity function rather than the Cole-Cole function is its numerical advantage in the convolutional algorithm. With the SE conductivity, we do not need to convert $\sigma(\omega)$ within each discretized voxel to $\sigma(t)$ because the SE conductivity has an explicit form in the time domain. The SE conductivity will not be beneficial when Maxwell’s equations are solved in frequency-domain, and we believe that is the reason why the SE conductivity has not been used extensively, except for the latest simulation study from Belliveau & Haber (2018).

Although the SE conductivity is not exactly the same as the Cole-Cole conductivity (eq. 1), their time-features are very similar, and when $c_{cc}$=1 (Debye model), they are equivalent. To illustrate cases when $c_{cc}$ is not equal to 1, we fit the Cole-Cole conductivity with the SE conductivity in time-domain; here, we update all three SE parameters: $\eta_{se}$, $\tau_{se}$, $c_{se}$ to fit Cole-Cole conductivity. Fig. 1 shows example Cole-Cole conductivity decays (t $>$ 0) with variable $c_{cc}$, and their fits with the SE conductivity. For the range of times of interest (10^-3-10¹ms), the SE function effectively fits the Cole-Cole, as shown in Fig. 1. They are essentially coincident. The estimated values of $\eta_{se}$ and $c_{se}$ are slightly smaller than their respective Cole-Cole counterparts; $\tau_{se}$ is coincident with $\tau_{cc}$ except when $c_{cc}$=0.2 (Table 1). Therefore, when interpreting the SE parameters, readers can use their understanding of Cole-Cole parameters and treat the SE and CC parameters as being similar. Note that we have used the impulse response of the Cole-Cole and SE functions when generating the fits. There is no loss of generality in doing this since the response due to an arbitrary waveform can be represented as a linear combination of impulse responses.

3Numerical experiments¶

To answer the four questions posed previously, we carry out TEM simulations using the SimPEG-EMIP code. For the spatial discretization, we use the 2D cylindrically symmetric mesh because of the cylindrical symmetry in the time-domain AEM system, which uses a horizontal loop (See Fig. 2). Rather than using a waveform for a specific AEM system, we use a step-off waveform as an input current. Again, there is no loss of generality since the response from an arbitrary waveform can be generated by a linear combination of step-off responses Fitterman & Anderson (1987). A horizontal receiver loop measuring the voltage (equivalent to -$db_z/dt$) is coincident with the source loop. A chargeable cylinder is embedded in the resistive halfspace ($\sigma_{half}$=10^-3 S/m). The depth to the top ($z_{top}$), radius ($r$), and thickness ($h$) of the chargeable cylinder are correspondingly $z_{top}$=50 m, $r$=200 m, and $h$=100 m; the SE parameters of the cylinder are $\siginf$=0.1 S/m, $\eta_{se}$=0.1, $\tau_{se}$=1ms, $c_{se}$=0.7; the cylinder is 100 times more conductive than the halfspace and its effective conductance ($\sigma h$) is 10 S.

To understand how negative transients are caused by the presence of chargeable rocks, we first explore how the electric field diffuses into the earth after the input current is turned off. Fig. 3(a) shows the simulated electric field in y-direction (into the page) at four different time channels (0.01-50ms). At early times (0.01-0.3ms) electric fields, which rotate in the horizontal plane in a counter-clockwise direction, are induced in both the halfspace and in the conductive cylinder. As time passes, the electric field diffuses downwards and radially outwards; particularly large rotating electric fields are induced in the conductor. These inductive currents are responsible for “charging up” the earth material. At a later time (6 ms), the inductive currents have gone and only the decaying polarization currents remain. The resultant electric fields (and currents) have reversed direction; this is due to IP effects. At a later time (50ms), these IP effects have decayed away. With Faraday’s law, electric fields are generated by time-varying magnetic field ($d\b/dt$), and the measured voltage is the same as the vertical component of -$d\b/dt$. Similarly, Fig. 3(b) shows the vectoral distribution of $d\b/dt$ in time. At 0.3ms, the high amplitude of $d\b/dt$ is shown in the target, and the main direction of $d\b/dt$ (white arrow) is downward. However, at 6 ms the upward direction of $d\b/dt$ (red arrows) is generated by IP effects; this will result in negative transients at the receiver loop. It is important to notice that electric fields generated either from EM or IP effects (Fig. 3a) do not cross a boundary. There is no charge build up on the boundary and there are no channeled currents (which is the mechanism by which IP signals for the grounded DC-IP surveys are generated). The current channeling could happen if the cylindrical symmetry is broken (e.g. the source loop is located away from the center of the chargeable cylinder), but our analysis is focused on when cylindrical symmetry is preserved; the IP effects we show are solely due to the inductive polarized currents.

Based upon the physical understanding of IP effects due to a loop source, we examine the data measured at the receiver loop. The black lines in Fig. 4(a) show the measured time decays, $d^{obs}$, (on a log-log scale); negative values are shown after 2 ms (black dashed line). Another simulation is carried out without IP effects ($\eta_{se}$=0), and the computed data are shown with the blue line (no negatives); we call these the fundamental data, $d^F$; they include only EM induction effects. The IP data, $\dip$, are defined as

The system noise-level is set to 10^-4 pV/A-m² based upon a field data set measured at Mt. Milligan with a VTEM system (Figure 4.23 in Kang (2018)), which denoted as the grey shaded region in Fig. 4(a). At the early times ($t<$1 ms), $d^{obs}$ and $d^{F}$ are almost coincident indicating that EM induction dominates the response. On the other hand, IP effects are dominant at later times ($t>$ 2ms). To show the relative strength of the IP effects, we define the ratio, $R$, between $|d^F|$ and $|\dip|$:

We show $R$ in Fig. 4(b); ratios smaller than 10^-2 are ignored. Between 1 ms and 40ms $R$ is greater than 0.1, indicating that there are considerable IP effects in the observations. When $R$=1, the observation is zero which corresponds to the time that the sign reversal occurs.

The EM induction processes within the background conductivity structure influence the electric field, which serves as a forcing function for IP effects. These IP effects translate into the reversed direction of the electric field in the target which results in negative transients that are observed after the EM induction effects have decayed. In the following sections, we carry out TEM-IP simulations to systematically investigate the feasibility of the airborne IP technique. Variable model parameters are shown in Fig. 2. Considering the typical AEM system specifications, we limit our attention to the measured time range from 10^-2ms to 10¹ms and to voltages greater than the noise level (10^-4 pV/A-m²). Further, for the metric pertaining to whether we can see IP signals or not, we use the existence of the negative datum being greater than the noise floor in the measured time range. Namely, we ignore subtle IP signals smaller than EM signals ($R<1$).

3.1Effects of background conductivity on IP signals¶

As we demonstrated in Fig. 3, the electric fields are the forcing functions which generate IP signals. Understanding the effects of the background conductivity on these electric fields is therefore crucial for understanding the resultant IP signals. To investigate this, we perform TEM-IP simulations for a range of values of $\siginf$ and $\sigma_{half}$. Other parameters for the simulation setup are the same as those in Fig. 1 ($r$=200m, $h$=100m, $z_{top}$=50m, $\eta_{se}$=0.1, $\tau_{se}$=1ms, $c_{se}$=0.7).

The first experiment involves varying $\siginf$, which ranges from 10^-4 S/m to 1 S/m, while $\sigma_{half}$ is fixed to 10^-3 S/m. Figs. 5(a) and (b) show the simulated time decays and corresponding ratios, $R$. Negative transients are only visible when $\siginf$ is 0.01 S/m and 0.1 S/m (red and green curves). When $\siginf$ is too high (e.g. 1 S/m), EM effects dominate at all times and no negative transients are visible in the observed data. At the other end of the spectrum, the very resistive target shows the smallest $R$. These results show that a moderately conductive target provides the best opportunity for observing negative transients in the data. When $\siginf$=1 S/m, there are no negatives in the time decay curve.

To explore the effect of the halfspace conductivity, $\sigma_{half}$, we fix the ratio $\siginf$/ $\sigma_{half}$ to be 10, and change $\sigma_{half}$ from 0.1 S/m to 10^-4 S/m. In Fig. 6, negatives are present when $\sigma_{half}$ is 10^-4 S/m and 10^-3 S/m, but not for the other cases. This shows that when the conductivity of the non-chargeable halfspace is too high ($>$0.1 S/m), measuring IP signals will be challenging even though chargeable materials exist. Therefore, a moderately conductive target ($\sim$0.01-0.1 S/m) in a resistive host ($\sim$10^-4-10^-3 S/m) provides the best circumstances for observing strong IP signals.

3.2To what depth can we expect to detect chargeable material with airborne IP?¶

Often, the maximum depth that airborne IP can see chargeable targets is considered to be fairly low ($z<\sim$100 m from Macnae (2016)). However, the possibility exists to see deeper when the host is resistive and the chargeable target is moderately conductive. Here we explore detectability of a chargeable target by altering the depth of the target ($z_{top}$) from 0m to 350 m. Fig. 7 shows the time decays with variable $z_{top}$ when $\siginf$ and $\sigma_{half}$ are 0.1 S/m and 10^-3 S/m respectively. Negatives are present when $z_{top}$ ≤ 200m. By decreasing $\sigma_{half}$, this depth can be increased to 300 m as shown in Fig. 8. Hence, it is possible to detect a deeper chargeable target using the airborne IP technique when the target is moderately conductive and the host rock is resistive (10^-4 S/m). For instance, at the Tli Kwi Cho kimberlite deposit, negatives were measured near a kimberlite pipe. This moderately conductive pipe was embedded in a resistive host rock (10^-4 S/m); it was located $\sim$70 m below the surface and its radius and thickness were approximately 150 m and 200 m, respectively Kang et al. (2017). This geometry is similar to that of our chargeable cylinder shown in Fig. 2.

The maximum depth that we can see negatives will depend upon IP parameters. For instance, greater chargeability will increase the strength of the IP signals and therefore the maximum depth can be increased with increased chargeability Macnae (2016). The effects of the time constant are more complicated to unravel. We explore this by changing the time constant of the target from 1ms to 10ms and altering the depth of burial. We obtain Fig. 9 and observe that the maximum depth is decreased from 350m to 250m. Performing a similar analyses for the time constant ranging from 0.1ms to 10s, we obtain the maximum depth as a function of the time constant as shown in Fig. 10. The maximum depth starts from zero when $\tau_{se}$=0.1ms, increases until $\tau_{se}$=1ms, and then decreases as $\tau_{se}$ increases. We simulated two cases in which $\sigma_{half}$ was 10^-3 S/m and 10^-4 S/m but the conductivity of the target ($\siginf$) was fixed to 10^-1 S/m. Greater maximum depth is shown when $\sigma_{half}$=10^-4 S/m. Hence there is an optimal time constant ($\sim$1ms) that can generate the greatest IP signals. This can be understood from the following. Considering the measured time range of the data: 10^-2-10ms, there simply not enough time to charge up material that has time constant greater than 3s. Further, when the IP decay is too fast (small time constant) compared to EM decay, it is hard to be the signal in the observation.

3.3The effects of target size¶

The examples shown so far have illustrated general principles concerned with the ability to detect IP bodies at depth. We have dealt with a specific geometry and have worked with a fairly large target. The strength of the IP signal depends upon the size and geometry of the target. To begin an exploration of the impact of target size on detectability, we show the effects of making the body smaller. We first reduce the radius ($r$) from 200 m to 50 m. As a result, the maximum depth has decreased from 300m (Fig. 8) to 100 m as shown in Fig. 11. We can also reduce the thickness ($h$), as shown in Fig. 12; here $z_{top}$=150m and $r$=100m. When $h$=10m, we no longer observe the negatives. Depending on the geological setting, more complicated situations may occur and the potential for seeing an IP signal in the airborne data will require 3D modelling appropriate to the geology. For instance, layering of the subsurface can also make significant impact to the maximum depth in practice, and this was not taken account in our analyses.

3.4Which are the characteristics of detectable chargeable materials in AEM data?¶

The different polarization characteristics of earth materials can be translated into different Cole-Cole or SE parameters. In particular, their frequency spectrum, or the time period over which the polarization process is occurring, differs; this time period is closely related to the time constant, τ. For instance, fine-grained sulphides will show a higher frequency spectrum and smaller time constant than coarse-grained sulphides. Ground and airborne surveys such as DC-IP and AEM often use different base frequencies: e.g. 0.125 Hz and 25 Hz. As AEM surveys have a higher base frequency, they are more sensitive to chargeable materials which are characterized by high frequency or smaller time constants. Three main chargeable targets of interest in airborne IP are (a) fine-grained sulphides Pelton et al. (1978)Revil et al. (2017), (b) clays Macnae (2016)Leroy & Revil (2009), and (c) ice Grimm & Stillman (2015)Kang et al. (2017). In Fig. 13, we have plotted the frequency spectrum of each of these materials along with the frequency spectrum of DC-IP and AEM surveys.

To examine how each of these materials impacts the observed time-decays, we have defined four models in Table 2 for which we will simulate AEM data. The decays, plotted in Fig. 14, show the characteristic time behavior associated each of the chargeable materials:

• Type A: Typical time decay showing positive early time data and negative late time data; this can be generated by fine-grained sulphides and clays.

• Type B: Double sign-reversal; when sulphides or clays have very fine grain size the resulting time constant is smaller.

• Type C: No negatives, but a positive ‘bump’ at late time; this is when there is a deep conductor below a chargeable target, which can generate strong positive EM signals at late time.

• Type D: No positives; this can be generated by an extremely chargeable target such as ice ($\simeq$0.9) located very near surface, or not measuring early enough time channels.

4Conclusions¶

We have developed a convolutional TEM simulation code that directly solves Maxwell’s equations in time with the SE (stretched exponential) conductivity function. The SE conductivity is a good representation of the Cole-Cole conductivity for the typical time range used in AEM. With our simulations, we showed that:

• Negative transients in AEM data can be caused by the reversed direction of the electric field in a chargeable target and are visible at late times when EM induction is small.

• Moderately conductive targets (0.01-0.1 S/m) in a resistive host (10^-4-10^-3 S/m) show the best potential for generating strong IP signals (negatives) in AEM data.

• The depth at which we can detect a target with airborne IP depends upon the background conductivity, but for an ideal situation (a conductive, chargeable target in a resistive host), the target can be detected up to 300 m depth provided its time constant is close to 1ms. This maximum depth will naturally be affected by the layering, but this was not taken account in our analyses.

• The strength of the IP signals and the depth of detectability of a target is dependent upon the size of the target and its geometry. In general 3D simulations, cast within the relevant geologic context, are required.

• The three main sources of chargeable material detectable in AEM (fine grain sulphides, clay and ice) can give rise to four different characteristic decay curves.

The overriding question of practical concern is whether, or under what circumstances, you can see an IP target at depth in airborne EM data. The situation is complex and cannot be answered by a simple, fixed depth of investigation rule. Forward simulation which emulates the potential geology and its associated physical properties, is required. To advance this capability, we have developed the SimPEG-EMIP codes as a part of the open-source software project, SimPEG (https://www.simpeg.xyz). The main workhorse for this paper has been the 2D cylindrical mesh, and the figures in this text can be reproduced with provided Jupyter Notebooks (https://github.com/simpeg-research/kang-2018-AEM). Moreover, the source code is downloadable and users can explore the use of different parameters. We hope these efforts contribute to the wider challenge of extracting IP information from airborne time-domain EM data.