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

Abstract

The potential for extracting and interpreting induced polarization (IP) data from airborne surveys is now broadly recognized. There is, however, still considerable discussion about the conditions under which the technique can provide knowledge about the subsurface and thus, its practical applications. Foremost among these is whether, or under what conditions, airborne IP can detect chargeable bodies at depth. To investigate, we focus on data obtained from a coincident-loop time-domain system. Our analysis is expedited by using a stretched exponential rather than a Cole-Cole model to represent the IP phenomenon. Our paper begins with an example that illuminates the physical understanding about how negative transients (the typical signature of an IP signal in airborne data) can be generated. The effects of the background conductivity are investigated; this study shows that a moderately conductive and chargeable target in a resistive host is an ideal scenario for generating strong IP signals. We then examine the important topic of estimating the maximum depth of the chargeable target that can generate negative transients. Lastly, some common chargeable earth-materials are discussed and their typical IP time-domain features are analyzed. The results presented in this paper can be reproduced and further explored by accessing the provided Jupyter notebooks.

Keywords:Airborne ElectromagneticsInduced PolarisationElectromagnetic Geophysics

...

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):

σcc(ω)=σηccσ1+(ıωτcc)ccc\sigma_{cc}(\omega) = \siginf - \frac{\eta_{cc} \siginf}{1+(\imath \omega \tau_{cc})^{c_{cc}}}

where σ\siginf is the conductivity at infinite frequency, ηcc\eta_{cc} is the chargeability, τcc\tau_{cc} is the time constant (s), and cccc_{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<z< 100 m). Viezzoli et al. (2017), showed the potential that a deeper chargeable target, such as a mineral deposit (z>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\J, in the frequency domain, can be written as:

J=σ(ω)E\J = \sigma(\omega)\E

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

j=σ(t)e\j = \sigma(t) \otimes \e

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

×e=bt\curl \e = -\frac{\partial \b}{\partial t}
×μ1bj=js\curl \mu^{-1}\b - \j = \j_s

where b\b is the magnetic flux density (Wb/m2) and js\j_s (A/m2) 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, σ(t)\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, 1ustep(t)1-u_{step}(t), can be written as

σse(1ustep)=σ0(1ustep(t))σηseexp((tτse)cse)ustep(t)\sigma_{se} \otimes (1-u_{step}) = \sigma_0 \Big(1-u_{step}(t)\Big) -\siginf \eta_{se}\rm{exp}\Big(-\left(\frac{t}{\tau_{se}}\right)^{c_{se}}\Big)u_{step}(t)

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

t(σse(1ustep))=σse(t)δ(t)=σse(t)- \frac{\partial}{\partial t}\Big(\sigma_{se} \otimes (1-u_{step})\Big) = \sigma_{se}(t) \otimes \delta(t) = \sigma_{se}(t)

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

σse(t)=σδ(t)σηset1(tτse)cseexp((tτse)cse)ustep(t),\sigma_{se}(t)=\siginf \delta(t)-\siginf \eta_{se}t^{-1}\left(\frac{t}{\tau_{se}}\right)^{c_{se}}\rm{exp}\Big(-\left(\frac{t}{\tau_{se}}\right)^{c_{se}}\Big)u_{step}(t),

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 σ(t)\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 cccc_{cc}=1 (Debye model), they are equivalent. To illustrate cases when cccc_{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: ηse\eta_{se}, τse\tau_{se}, csec_{se} to fit Cole-Cole conductivity. Fig. 1 shows example Cole-Cole conductivity decays (t >> 0) with variable cccc_{cc}, and their fits with the SE conductivity. For the range of times of interest (10-3-101ms), the SE function effectively fits the Cole-Cole, as shown in Fig. 1. They are essentially coincident. The estimated values of ηse\eta_{se} and csec_{se} are slightly smaller than their respective Cole-Cole counterparts; τse\tau_{se} is coincident with τcc\tau_{cc} except when cccc_{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.

The stretched exponential (SE) fit of the Cole-Cole conductivity in the time domain. Solid lines are the impulse response using a Cole-Cole representation and the circles denote the SE response. Black, blue, and red colors correspondingly indicate when c_{cc} is 0.3, 0.5, and 0.7.

Figure 1:The stretched exponential (SE) fit of the Cole-Cole conductivity in the time domain. Solid lines are the impulse response using a Cole-Cole representation and the circles denote the SE response. Black, blue, and red colors correspondingly indicate when cccc_{cc} is 0.3, 0.5, and 0.7.

Table 1:Comparison of the Cole-Cole (CC) and the resultant Stretched Exponential (SE) parameters for variable cccc_{cc}. The chargeability, time constant, and the frequency component are correspondingly represented as η, τ, cc; subscripts CC and SE denote the Cole-Cole and SE parameters.

CCSECCSECCSE
ηcc\eta_{cc}0.10.090.10.090.10.09
τcc\tau_{cc} (ms)10.81111
cccc_{cc}0.30.20.50.40.70.6

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 -dbz/dtdb_z/dt) is coincident with the source loop. A chargeable cylinder is embedded in the resistive halfspace (σhalf\sigma_{half}=10-3 S/m). The depth to the top (ztopz_{top}), radius (rr), and thickness (hh) of the chargeable cylinder are correspondingly ztopz_{top}=50 m, rr=200 m, and hh=100 m; the SE parameters of the cylinder are σ\siginf=0.1 S/m, ηse\eta_{se}=0.1, τse\tau_{se}=1ms, csec_{se}=0.7; the cylinder is 100 times more conductive than the halfspace and its effective conductance (σh\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 (db/dtd\b/dt), and the measured voltage is the same as the vertical component of -db/dtd\b/dt. Similarly, Fig. 3(b) shows the vectoral distribution of db/dtd\b/dt in time. At 0.3ms, the high amplitude of db/dtd\b/dt is shown in the target, and the main direction of db/dtd\b/dt (white arrow) is downward. However, at 6 ms the upward direction of db/dtd\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, dobsd^{obs}, (on a log-log scale); negative values are shown after 2 ms (black dashed line). Another simulation is carried out without IP effects (ηse\eta_{se}=0), and the computed data are shown with the blue line (no negatives); we call these the fundamental data, dFd^F; they include only EM induction effects. The IP data, dIP\dip, are defined as

dIP=dobsdF.\dip = d^{obs} - d^F.

The system noise-level is set to 10-4 pV/A-m2 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<t<1 ms), dobsd^{obs} and dFd^{F} are almost coincident indicating that EM induction dominates the response. On the other hand, IP effects are dominant at later times (t>t> 2ms). To show the relative strength of the IP effects, we define the ratio, RR, between dF|d^F| and dIP|\dip|:

R=dIPdF.R = \frac{|\dip|}{|d^F|}.

We show RR in Fig. 4(b); ratios smaller than 10-2 are ignored. Between 1 ms and 40ms RR is greater than 0.1, indicating that there are considerable IP effects in the observations. When RR=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 101ms and to voltages greater than the noise level (10-4 pV/A-m2). 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<1R<1).

A chargeable cylinder embedded in a halfspace earth. The 13m-radius source loop is located 30 m above the surface. The depth to the top of the prism is denoted by z_{top}. Right: the 2D cylindrically symmetric mesh is used for TEM simulations.

Figure 2:A chargeable cylinder embedded in a halfspace earth. The 13m-radius source loop is located 30 m above the surface. The depth to the top of the prism is denoted by ztopz_{top}. Right: the 2D cylindrically symmetric mesh is used for TEM simulations.

EM fields diffusing in the earth: (a) electric field in y-direction and (b) time derivative of the magnetic field (d\b/dt); the chargeable cylinder is outlined in white. At early times (0.01-0.3 ms), EM induction is dominant; this results in inductive electric fields rotating counter-clockwise and d\b/dt fields going upward. However, at 6ms the direction of the electric field is reversed (clockwise) as a result of the chargeable cylinder; similarly d\b/dt fields go upward (red arrows); This results in negative transients at the receiver loop (See Fig. %s).

Figure 3:EM fields diffusing in the earth: (a) electric field in yy-direction and (b) time derivative of the magnetic field (db/dtd\b/dt); the chargeable cylinder is outlined in white. At early times (0.01-0.3 ms), EM induction is dominant; this results in inductive electric fields rotating counter-clockwise and db/dtd\b/dt fields going upward. However, at 6ms the direction of the electric field is reversed (clockwise) as a result of the chargeable cylinder; similarly db/dtd\b/dt fields go upward (red arrows); This results in negative transients at the receiver loop (See Fig. 4).

(a) Simulated TEM responses; d^{obs} (black) are observations (EM+IP), d^F (blue) are fundamental (EM) data, and \dip (red) are IP data. Solid and dashed lines distinguish positive and negative values. Signals beneath the noise level (10^{-4} pV/A-m^2) are shown within the grey zone. (b) R =\frac{|\dip|}{|d^F|} shows the relative strength of the IP signals compared to the fundamental induction effects. In the grey region R<1 and the strength of the IP signal is smaller than the EM signal.

Figure 4:(a) Simulated TEM responses; dobsd^{obs} (black) are observations (EM+IP), dFd^F (blue) are fundamental (EM) data, and dIP\dip (red) are IP data. Solid and dashed lines distinguish positive and negative values. Signals beneath the noise level (10-4 pV/A-m2) are shown within the grey zone. (b) R=dIPdFR =\frac{|\dip|}{|d^F|} shows the relative strength of the IP signals compared to the fundamental induction effects. In the grey region R<1R<1 and the strength of the IP signal is smaller than the EM signal.

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 σhalf\sigma_{half}. Other parameters for the simulation setup are the same as those in Fig. 1 (rr=200m, hh=100m, ztopz_{top}=50m, ηse\eta_{se}=0.1, τse\tau_{se}=1ms, csec_{se}=0.7).

The first experiment involves varying σ\siginf, which ranges from 10-4 S/m to 1 S/m, while σhalf\sigma_{half} is fixed to 10-3 S/m. Figs. 5(a) and (b) show the simulated time decays and corresponding ratios, RR. 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 RR. 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, σhalf\sigma_{half}, we fix the ratio σ\siginf/ σhalf\sigma_{half} to be 10, and change σhalf\sigma_{half} from 0.1 S/m to 10-4 S/m. In Fig. 6, negatives are present when σhalf\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 (\sim0.01-0.1 S/m) in a resistive host (\sim10-4-10-3 S/m) provides the best circumstances for observing strong IP signals.

(a) Time decay curves with a variable conductivity of the chargeable cylinder. Solid and dashed lines distinguish positive and negative values. The halfspace conductivity is fixed at 10^{-3} S/m, whereas \siginf varies (10^{-4}-0.1 S/m). (b) Plots of R=|\dip|/|d^F| . In the grey region R<1 and the strength of the IP signal is smaller than the EM signal. The legend for both plots is shown in (b).

Figure 5:(a) Time decay curves with a variable conductivity of the chargeable cylinder. Solid and dashed lines distinguish positive and negative values. The halfspace conductivity is fixed at 10-3 S/m, whereas σ\siginf varies (10-4-0.1 S/m). (b) Plots of R=dIP/dFR=|\dip|/|d^F| . In the grey region R<1R<1 and the strength of the IP signal is smaller than the EM signal. The legend for both plots is shown in (b).

(a) Time decay curves with a variable conductivity of the halfspace, \sigma_{half}. Solid and dashed lines distinguish positive and negative values. The ratio between the halfspace conductivity and the target conductivity ( \siginf/\sigma_{half}), is fixed to 10, whereas \sigma_{half} varies (10^{-4}-0.1 S/m). (b) Plots of the |\dip|/|d^F| ratio, R. In the grey region R<1 and the strength of the IP signal is smaller than the EM signal. The legend for both plots is shown in (b).

Figure 6:(a) Time decay curves with a variable conductivity of the halfspace, σhalf\sigma_{half}. Solid and dashed lines distinguish positive and negative values. The ratio between the halfspace conductivity and the target conductivity ( σ\siginf/σhalf\sigma_{half}), is fixed to 10, whereas σhalf\sigma_{half} varies (10-4-0.1 S/m). (b) Plots of the dIP/dF|\dip|/|d^F| ratio, RR. In the grey region R<1R<1 and the strength of the IP signal is smaller than the EM signal. The legend for both plots is shown in (b).

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<z<\sim100 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 (ztopz_{top}) from 0m to 350 m. Fig. 7 shows the time decays with variable ztopz_{top} when σ\siginf and σhalf\sigma_{half} are 0.1 S/m and 10-3 S/m respectively. Negatives are present when ztopz_{top} 200m. By decreasing σhalf\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 \sim70 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 τse\tau_{se}=0.1ms, increases until τse\tau_{se}=1ms, and then decreases as τse\tau_{se} increases. We simulated two cases in which σhalf\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 σhalf\sigma_{half}=10-4 S/m. Hence there is an optimal time constant (\sim1ms) 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.

(a) Time decay curves with a variable target depth (z_{top}) ranging from 0-350 m. The halfspace conductivity is 10^{-3} S/m. (b) Plots of R=|\dip|/|d^F| . In the grey region R<1 and the strength of the IP signal is smaller than the EM signal. The legend for both plots is shown in (b).

Figure 7:(a) Time decay curves with a variable target depth (ztopz_{top}) ranging from 0-350 m. The halfspace conductivity is 10-3 S/m. (b) Plots of R=dIP/dFR=|\dip|/|d^F| . In the grey region R<1R<1 and the strength of the IP signal is smaller than the EM signal. The legend for both plots is shown in (b).

(a) Time decay curves with a variable target depth (z_{top}) ranging from 0-350 m. The halfspace conductivity is decreased to 10^{-4} S/m compared to Fig. %s. (b) Plots of R=|\dip|/|d^F| . In the grey region R<1 and the strength of the IP signal is smaller than the EM signal. The legend for both plots is shown in (b).

Figure 8:(a) Time decay curves with a variable target depth (ztopz_{top}) ranging from 0-350 m. The halfspace conductivity is decreased to 10-4 S/m compared to Fig. 7. (b) Plots of R=dIP/dFR=|\dip|/|d^F| . In the grey region R<1R<1 and the strength of the IP signal is smaller than the EM signal. The legend for both plots is shown in (b).

(a) Time decay curves with a variable target depth (z_{top}) ranging from 0-350 m. The time constant (\tau) is increased to 10^{-2}s. (b) |\dip|/|d^F| ratio, R. Grey region indicates R<1 meaning the strength of the IP signal is smaller than EM signal.

Figure 9:(a) Time decay curves with a variable target depth (ztopz_{top}) ranging from 0-350 m. The time constant (τ) is increased to 10-2s. (b) dIP/dF|\dip|/|d^F| ratio, RR. Grey region indicates R<1R<1 meaning the strength of the IP signal is smaller than EM signal.

Maximum depth of the chargeable cylinder at which we can observe negative transients as a function of the time constant (ms). Light blue and orange lines indicate two conductivity models having different halfspace conductivity (\sigma_{half}): 10^{-4} S/m and 10^{-3} S/m, respectively.

Figure 10:Maximum depth of the chargeable cylinder at which we can observe negative transients as a function of the time constant (ms). Light blue and orange lines indicate two conductivity models having different halfspace conductivity (σhalf\sigma_{half}): 10-4 S/m and 10-3 S/m, respectively.

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 (rr) 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 (hh), as shown in Fig. 12; here ztopz_{top}=150m and rr=100m. When hh=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.

(a) Time decay curves with a variable target depth (z_{top}). The radius (r) of the chargeable cylinder is decreased from 200 m to 100 m compared to Fig. 8. (b) |\dip|/|d^F| ratio, R. Grey region indicates R<1 meaning the strength of the IP signal is smaller than EM signal.

Figure 11:(a) Time decay curves with a variable target depth (ztopz_{top}). The radius (rr) of the chargeable cylinder is decreased from 200 m to 100 m compared to Fig. 8. (b) dIP/dF|\dip|/|d^F| ratio, RR. Grey region indicates R<1R<1 meaning the strength of the IP signal is smaller than EM signal.

(a) Time decay curves with a variable target thickness (h). The radius (r) and depth to the top (z_{top}) of the chargeable cylinder are fixed to 50 m and 150 m, respectively, whereas h varies (10-100 m) as shown in the legend in (b) which displays the |\dip|/|d^F| ratio. Grey region indicates R<1 meaning the strength of the IP signal is smaller than EM signal.

Figure 12:(a) Time decay curves with a variable target thickness (h). The radius (rr) and depth to the top (ztopz_{top}) of the chargeable cylinder are fixed to 50 m and 150 m, respectively, whereas hh varies (10-100 m) as shown in the legend in (b) which displays the dIP/dF|\dip|/|d^F| ratio. Grey region indicates R<1R<1 meaning the strength of the IP signal is smaller than EM signal.

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 (\simeq0.9) located very near surface, or not measuring early enough time channels.

Frequency spectrum of EM systems (DC-IP and AEM) and various IP sources (sulphide, clays, and ice). The time constant is denoted at the top x-axis. This figure is based upon previous research

Figure 13:Frequency spectrum of EM systems (DC-IP and AEM) and various IP sources (sulphide, clays, and ice). The time constant is denoted at the top xx-axis. This figure is based upon previous research Pelton et al. (1978)Revil (2013)Grimm & Stillman (2015)Macnae & Hine (2016)

.

Four different types of time decays (A-D) from different sources of IP. Parameters used to compute time decays are summarized in Table %s.

Figure 14:Four different types of time decays (A-D) from different sources of IP. Parameters used to compute time decays are summarized in Table 2.

Table 2:Parameters of the chargeable cylinder for the type curves (A-D) shown in Fig. 9. For Type C a conductive layer (0.1 S/m) is added 300 m below the surface; its thickness is 100 m.

DivisionType AType BType CType D
LithologyClay
Sulphide
Clay (finer)
Sulphide (finer)
Type A
with a deep conductor
Ice
ztopz_{top} (m)5050500
σhalf\sigma_{half} (S/m)10-310-410-310-4
σ\siginf (S/m)2×10-22×10-22×10-210-3
ηse\eta_{se}0.10.10.10.9
τse(ms)\tau_{se} (ms)10.110.08
csec_{se}0.70.70.70.5

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.

References
  1. Cole, K. S., & Cole, R. H. (1941). Dispersion and Absorption in Dielectrics I. Alternating Current Characteristics. The Journal of Chemical Physics, 9(4).
  2. Pelton, W., Ward, S., Hallof, P., Sill, W., & Nelson, P. (1978). Mineral discriminatino and removal of inductive coupling with multifrequency IP. Geophysics, 43(3), 588–609.
  3. Tarasov, A., & Titov, K. (2013). On the use of the Cole-Cole equations in spectral induced: Polarization. I, 195(1), 352–356.
  4. Fink, J., McAlister, E., Sternberg, B., Wieduwilt, W., & Ward, S. (1990). Induced Polarization Applications and Case Histories (J. B. Fink, E. O. McAlister, B. K. Sternberg, W. G. Wieduwilt, & S. H. Ward, Eds.). Society of Exploration Geophysicists.
  5. Grimm, R. E., & Stillman, D. E. (2015). Field Test of Detection and Characterization of Subsurface Ice Using Broadband Spectral Induced Polarization. 38(November 2014), 28–38.