Spectral Induced Polarization Survey with Distributed Array System for Mineral Exploration: Case Study in Saudi Arabia

Abstract

The Saudi Arabian Glass Earth Pilot Project is a geophysical exploration program to explore the upper crust of the Kingdom for minerals, groundwater, and geothermal resources as well as strictly academic investigations. The project began with over 8000 km² of green-field area. Airborne geophysics including electromagnetic (EM), magnetics, and gravity were used to develop several high priority targets for ground follow-up. Based on the results of airborne survey, a spectral induced polarization (SIP) survey was completed over one of the prospective targets. The field data were collected with a distributed array system, which has the potential for strong inductive coupling. This was examined in a synthetic study, and it was determined that with the geometries and conductivities in the field survey, the inductive coupling effect may be visible in the data. In this study, we also confirmed that time domain is vastly superior to frequency domain for avoiding inductive coupling, that measuring decays from 50 ms to 2 s allow discrimination of time constants from 1 ms to 5 s, and the relaxation parameter C is strongly coupled to intrinsic chargeability. We developed a method to fully include all 3D EM effects in the inversion of induced polarization (IP) data. The field SIP data were inverted using the generalized effective-medium theory of induced polarization (GEMTIP) in conjunction with an integral equation-based modeling and inversion methods. These methods can replicate all inductive coupling and EM effects, which removes one significant barrier to inversion of large bandwidth spectral IP data. The results of this inversion were interpreted and compared with results of drill hole set up in the survey area. The drill hole intersected significant mineralization which is currently being further investigated. The project can be considered a technical success, validating the methods and effective-medium inversion technique used for the project.

1. Introduction

The Saudi Arabian Glass Earth Pilot Project is a multiyear Earth observation program. The goal is to explore for resources in the relatively unexplored regions of Saudi Arabia for minerals, groundwater, and geothermal resources as well as strictly academic investigations. The initial stage of the project involved the acquisition, processing, and interpretation of airborne electromagnetic, gravity, and magnetic geophysical data over an 8000 square kilometer area [1]. The project helped to identify a dozen potential mining targets, which can be associated with gold and other base metal deposits, including zinc and copper. A follow-up ground spectral induced polarization (SIP) survey was proposed to further investigate one of the identified anomalies. The SIP data were collected by Dias Geophysical with their distributed array system.

The spectral induced polarization (SIP) data collected in 3D present a very rich data set. Modern distributed array systems are instruments which allow nearly arbitrary survey geometry. They can collect 10,000 s of sounding points over several decades of bandwidth points (e.g., [2,3]). These data sets contain information about the 3D distribution of time dependent or frequency dependent complex conductivity of the medium. The advantage of this information is that it can be translated into knowledge about mineralization which goes well beyond that contained in conventional conductivity. In fact, it may be possible to discriminate economic and non-economic mineralization from the surface using this technique (e.g., [4,5]).

However, the major challenges to a widespread use of this technology were a lack of a model which describes and properly characterizes the complex conductivity response of different media (microscopic effects), and a lack of suitable software and methods to process and interpret large volumes of the SIP data with a technique which is rigorous enough to honor all the macroscopic physical effects. We address these two issues with our new inversion method and code which use the new induced polarization model: the generalized effective medium theory of induced polarization (GEMTIP) introduced by Zhdanov [5]. The inversion code also honors all the macroscopic physics of electromagnetic field which is required for modeling.

We first consider the physical model describing the induced polarization phenomenon. The induced polarization is related to a complex (amplitude and phase or real and imaginary) electrical resistivity which varies with frequency. Many different models have been proposed to replicate the complex resistivity spectra seen in real rocks. These include many models such as Debye dispersion [6], but most prevalent is the Cole–Cole model [7]. It was introduced to geophysics by [8,9], and [10]. A review of this method was done by [11]. Indeed, the Cole–Cole model is the most widely applied model for SIP data (e.g., [12,13,14,15,16]). However, it does not connect the mineralization and rock composition directly to the spectral induced polarization spectrum.

In 2008, Zhdanov introduced the GEMTIP model which does allow direct connection of the rock matrix to the complex resistivity spectrum [5]. In this paper, Zhdanov showed that the Cole–Cole model could be derived analytically from the general effective medium theory for electromagnetic field. Indeed, the Cole–Cole model is a specific solution of the GEMTIP model when conductive inclusions are spherical and there is only one type or phase of a mineral which 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 [17,18], who used QEMScan mineralogical testing along with laboratory measurements of the complex spectra of the samples.

There are also challenges with modeling all EM induction and IP effects accurately and then inverting the data on a large scale. To simplify the problem, most inversion methods are based on a direct current (DC) approximation to Maxwell’s equations, and then a linear approximation of the IP effects with conductivity [19,20]. The problem with this approach is three-fold: (1) inductive coupling (IC) between the transmitter and receiver wires and also within the earth are ignored while these effects can be large [21]; (2) the assumption that IP effects are linearly related to conductivity is often invalid; and (3) other IP parameters (time constant, GEMTIP parameters, etc.) cannot be included using this method. Other researchers have included the full solution to Maxwell’s equations to eliminate the traditional limitations to inductive coupling. For example, [22] uses a full solution but simplified the solver for a single in-phase and quadrature resistivity. The authors of [23] approached this problem with a full Cole–Cole relaxation model and accurate transmitter waveform and inductive coupling modeling, but the method was of limited applicability because the solution was restricted to 1D conductivity structures. [24] used a half-space approximation to fit the IC and then remove this part before inversion. Others have modeled the full electromagnetic (EM) and IP response for a 3D conductive body for ground-based methods (e.g., [25]) and airborne EM (e.g., [26]). But no one has previously solved the full EM problem in 3D during inversion to include inductive coupling in a ground-based IP survey.

The full solution to this problem is difficult. Computational times and memory requirements are problematic, especially when including full wire paths with large 3D data sets and 3D earth models. This is exacerbated because many forward modeling runs are required for each inversion run. Also, because of the ill-posedness of the problem [27,28,29], many inversion runs are required to properly explore the range of models which are acceptable to the observed data.

These problems have been overcome in part by using a very efficient integral equation solution [30] and a re-weighted regularized conjugate method to keep memory requirements at a minimum. Two papers have been written on this inversion method which includes both full IC and determination of all four Cole–Cole parameters: [15,31]. This work was extended by [32] to apply full 3D inversion including EM effects and complex transmitter wire paths, and then to solve for four parameters of the GEMTIP model.

In this paper, we applied these advanced modeling and inversion methods to the spectral IP data collected in Saudi Arabia during the Glass Earth Pilot Project. We used an efficient integral equation forward modeling scheme, including full EM effects and wire paths, for the forward modeling. The earth is represented as voxels with the complex conductivity of each voxel being parameterized with the GEMTIP model. We solved for the four dominate GEMTIP parameters for each voxel. The non-uniqueness of the problem is mitigated with regularized inversion. The re-weighted conjugate gradient method was used to optimize the search for the minimum of the objective functional.

Based on the results of the inversion, the potential target was identified in a green-field area of Saudi Arabia. The target was drilled, and mineralization was confirmed in a drill hole representing a significant technical success.

GEMTIP Resistivity Relaxation Model

For completeness, we begin our discussion with a brief review of the GEMTIP resistivity relaxation model. In a general case, the effective conductivity of rocks is not necessarily a constant and real number but is complex and may vary with frequency. A general approach to constructing the resistivity relaxation model is based on the rock physics and description of the medium as a composite heterogeneous multiphase formation. The complex resistivity is parameterized using the GEMTIP model [5], which is based on effective-medium representation of a composite heterogeneous multi-phase rock formation. The GEMTIP model was developed to characterize the complex resistivity of multiphase heterogeneous rocks and their petrophysical and structural properties, including grain size, grain shape, porosity, anisotropy, polarizability, volume fraction, and conductivity of the inclusions in the pore space.

In the papers [17,18], we have demonstrated that the GEMTIP model of the complex resistivity for a three-phase medium with elliptical inclusions can be described by the following formula:

$$\sigma \left(\omega \right)={\sigma }_o{\left(1+{\displaystyle \sum }_{l=1}^2{\displaystyle \sum }_{\alpha =x,y,z}^{}\frac{f_l}{3{\gamma }_{l\alpha }}\left[1-\frac{1}{1+s_{l\alpha }{\left(i\omega {\tau }_l\right)}^{C_l}}\right]\right)}^{}$$

(1)

where ${\sigma }_o$ is the matrix conductivity.

Parameters ${\tau }_l$ and $C_l$ are similar to the Cole–Cole model and represent the time constants and relaxation parameters, respectively. The coefficients ${\gamma }_{l\alpha }$ and $s_{l\alpha }$ are the structural parameters defined by the geometrical characteristics of the ellipsoidal inclusions, and $f_l$ is their volume fraction. In a special case of only spherical grains with one mineral phase in a background rock matrix, Equation (1) reduces to a Equation (2) which is very similar to the widely used Cole–Cole formula for complex resistivity [5]:

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

(2)

where ${\sigma }_o$ is the DC conductivity (S/m); $f$ is volume fraction of mineral phase; $\omega$ is the angular frequency (rad/s), $\tau$ is the time parameter; and C is the relaxation parameter.

Equations (1) and (2) for complex and frequency dependent conductivity make it possible to consider simultaneously the electromagnetic induction and IP effects in observed EM data within both forward modeling and inversion stages of interpretation. Therefore, in the framework of this approach we do not need to apply any special technique of separating those effects, which is often done in interpretation of the IP data.

2. Spectral IP Survey Description

2.1. Geological Background

There are three distinct geologic regimes in the original study area (Figure 1 and Figure 2). The western most part contains Arabian shield with mountains up to 1 km high. The central part is formed by Cenozoic volcanic flows that are several hundred meters thick and underlain by the Arabian shield. Eruptions were recorded as recently as 1256 AD. The eastern part of the survey area is the focus of this study. It exhibits Arabian shield with sediment filled basins and hills hundreds of meters high. The Quaternary sediments in the areas are predominately alluvium, Aeolian sand, and slightly consolidated alluvium of old wadi deposits. There are many historic mines in the areas, including the largest mine in Saudi Arabia, the Mahd Ad Dahab (“Cradle of Gold”). A geologic map of the survey area is shown in Figure 1. The legend for this figure is shown in Figure 2.

2.2. Geophysical Background

The original regional study [1] area was covered by airborne EM, magnetics, and gravity surveys. From the airborne data, several large-scale regional features were discovered [1]. Many small-scale anomalies (hundreds of meters areal extent) were also seen based on the Glass Earth models of the physical properties of the subsurface. Note that, most of these anomalies are located on or very close to the large regional features that have been identified by geophysical data. Sulfides, such as pyrite, pyrrhotite, and chalcopyrite, are commonly associated with gold and copper mineralization. These minerals are electrically conductive and can often be detected by airborne surveys. We believe some of these anomalies may be sourced from sulfides. The fact that both copper and gold mineralization is known in the surrounding area is very encouraging.

An example of one of these anomalies is identified in the airborne EM data shown in Figure 3. The targets that we choose to follow-up with are all in exposed shield, all are conductive anomalies, and many have associated magnetic signatures. Each anomaly has been checked for geophysical significance and is not related to any cultural artifacts, neither present in the satellite images, nor marked by the aircrew.

We made one of the selected targets the highest priority target. The reasons for this were because the shape of the target suggested alteration or mineralization, and not saline groundwater or overburden, which are likely very common conductive features in this region. Also, there are magnetic anomalies associated with the conductive anomaly, which is another strong sign for a possible mineralization. The conductivity associated with this deposit, as inferred from the inversion of airborne data, is shown in Figure 3.

From these images, it appears the body exhibits both a magnetic response and an elevated electrical conductivity. The body is cigar shaped dipping steeply to the south–southeast. It is located near a large regional structure mapped with both the gravity and electromagnetic methods. A conductive overburden can be seen as a near surface target, but the prospective mineralization body is clearly hosted within the resistive bedrock.

The outstanding question is “Does this body contain mineralization?”. This cannot be wholly established by geophysics alone but determining if the body has an induced polarization response increases the certainty that this conductivity anomaly is due to mineralization, and it is possibly an economic mineral deposit. The spectral induced polarization method is a great compliment to airborne EM. Being a resistive method, it does an excellent job of imaging targets of moderate resistivity or conductivity. It also can differentiate between disseminated sulfides, massive sulfides, and clays or groundwater of similar conductivity.

2.3. Field Data Collection

Two SIP surveys were collected over the same area. One was a conventional 2D dipole–dipole type setup with five lines spaced 350 m apart and 3 km long and 100 m dipoles, and another was a true 3D type survey. Dias Geophysical collected both spectral IP surveys with their 3D distributed array system [33] in February of 2018. The ground IP survey was designed with five three-kilometer-long lines oriented perpendicular to the dominant geologic strike in the area. The survey used 100 m grounded dipoles as sources and receivers to best collect spectral IP data. These lines are shown in Figure 4. The SIP data are collected as a time series of voltage with each current injection. The majority of the current injection dipoles used a 50% duty cycle 8 s square wave. Each current injection point undergoes several minutes of current injection with the waveform being recorded at the injection points. The voltage at the receivers is recorded continuously throughout their entire deployment.

The second survey was very similar to the first one, but had three lines of receivers, and a single line of transmitting dipoles along the center line of receivers. Voltages were measured at each receiver potential electrode for every transmitter position, allowing a large number of combinations of off-line and on-line data points to be measured for each transmitter position. This survey layout is also shown in Figure 4.

2.4. Field Data Processing

The data reduction stage converts the receiver voltages measured at a common ground to dipole voltages. The receiver voltage time series is then paired with the injection current, stacked to reduce noise, and then integrated through the time windows. A primary voltage measured during the current on-time is also extracted. The data were processed to 20 time windows from 40 ms to 1.92 s window centers, plus Vp (on-time). Before inversion, spikes, and outliers were removed. However, due the distributed array geometry, visual quality control of the data is difficult. For the wire paths which are nearly perpendicular to the transmitters at the endpoints, 3D features have a large effect on the response, and negative and positive apparent resistivities are possible. To deal with this, during inversion these responses were examined and the sign of apparent resistivities in these specific cases was corrected. Example maps and decays of the observed data are shown in the field data analysis section.

3. Analysis of Inductive Coupling (IC) Effect

3.1. Inductive Coupling in IP Measurements

The collection of true broadband spectral IP data requires higher frequencies or early times, which increases the chances of incurring coupling. Additionally, as surveys become more complex and truly 3D, it is more difficult to eliminate the IC effect. The inductive coupling response can be larger than the IP response and it has been a major hindrance in the past to interpretation. The IC effect has been examined several times before, mostly in the context of standard array configuration (e.g., [21]). It has also been attempted to remove the IC effects and to correct the data to pure IP and resistivity responses (e.g., [33]).

In the ground SIP survey design implemented in this project, only one common receiver wire needs to be used for all receiver locations, greatly simplifying survey logistics and increasing the data collection speed. Each receiver can be referenced to every other receiver in the processing to create many receiver dipoles. However, this produces long receiver wire paths, some of which run parallel to transmitter current dipoles, and hence inductive coupling can become problematic if not modeled correctly.

In this section, we look at the effects of inductive coupling and advanced survey designs in the context of extracting spectral IP information. In order to illustrate our approach, we first consider the distributed field survey layout shown in Figure 5. There are nine transmitter dipoles operated in the case, and each neighboring current electrode pair operated as a current bipole. The voltage at each receiver point is measured relative to a common ground.

A note on our nomenclature bipole and dipole: bipole is defined for our purposes as two poles connected by a length of wire. The separation between the poles is great enough that they cannot be considered at the same point. A current dipole is a theoretical construct in which the poles are proximal enough relative to the measurement distance that, for calculations, their separation distance approaches zero while the current strength approaches infinity in such a way to keep the moment constant.

We will demonstrate in this section that most of the IC effect is due to wire to wire, not EM effects in the earth. This means that simply using the full sets of Maxwell’s equations but not including the true wire path is not sufficient for accurate modeling. In fact, in most cases using the straight wire path is nearly equivalent to using a DC approximation.

3.2. Three-Dimensional Forward Modeling Using Integral Equation Method

The forward modeling for electric field is based on the integral equation method. This has been covered extensively in the literature (e.g., [28,29]). For our purposes we start with the expression for the total electric field at point j at the surface of the domain of interest:

$$\mathit{E}\left(r\prime ,\omega \right)=\underset{D}{\overset{}{{\displaystyle \iiint }}}{\mathit{G}}_E\left(r\prime ,r,\omega \right)\Delta \tilde{\sigma }\left(r,\omega \right)\mathit{E}\left(r,\omega \right)dV+{\mathit{E}}^b\left(r^{\prime },\omega \right)$$

(3)

where $r^{\prime }$ is an arbitrary point along the receiver wire; $\mathit{E}$ is the total electric field inside the domain of interest; ${\mathit{E}}^b$ is the background electric field produced by a 1D earth; ${\mathit{G}}_E$ is the domain-to-receiver electric Green’s tensor; and $\Delta \tilde{\sigma }$ is the anomalous conductivity, which is a departure from the 1D earth structure. The total electric field is a superposition of the anomalous field caused by the conductivity anomalies and the background field produced by the layered earth in the absence of any anomalous structure. Note the tilde over the anomalous conductivity indicating that this is a complex frequency dependent term to include the induced polarization effects, as described by Formula (1).

The electric field at the receiver’s midpoint is commonly used as the predicted data for modeling and inversion (e.g., [32]). This is typically a good approximation for dipole data, which only may differ from the actual data for the very near offsets (when distance between the transmitter and receiver dipoles is less than two times their length). In a dipole–dipole survey, the voltage is measured as a difference between two adjacent electrodes. The voltage can be accurately approximated as the inline directed electric field at the midpoint multiplied by the distance between the two points.

However, with modern distributed array systems, the electric field at a receiver midpoint multiplied by the dipole length may be a poor approximation to the actual measured voltage. For example, the voltage may be measured between two distinct points on the different lines (e.g., $M_2$ and $M_9$ as shown in Figure 5). It will be much more accurate to integrate the electric field along the wire paths to obtain the voltage as given by Equation (3) than to use a midpoint electric field approximation. The voltage for the $j\mathrm{th}$ receiver combination is thus synthesized for a measurement from electrode $M_k$ to electrode $M_l$ as follows:

$$\begin{gathered} \Delta V_j=\Delta V_{M_k{M}_l}\left(\omega \right)=\underset{M_k}{\overset{M_l}{{\displaystyle \int }}}\mathit{E}\left(r^{\prime },\omega \right)·\mathit{d}\mathit{l}= \\ \begin{array}{c}\underset{M_k}{\overset{M_l}{{\displaystyle \int }}}\left(\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\mathit{G}}_E\left(r^{\prime },r_i,\omega \right)\Delta \tilde{\sigma }\left(r,\omega \right)\mathit{E}\left(r_i,\omega \right)dV+{\mathit{E}}^b\left(r^{\prime },\omega \right)\right)·\mathit{d}\mathit{l}=\\ \underset{M_k}{\overset{M_l}{{\displaystyle \int }}}\left(\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\mathit{G}}_E\left(r^{\prime },r_i,\omega \right)\Delta \tilde{\sigma }\left(r_i,\omega \right)\mathit{E}\left(r_i,\omega \right)dV\right)·\mathit{d}\mathit{l}+\underset{M_k}{\overset{M_l}{{\displaystyle \int }}}{\mathit{E}}^{\mathit{b}}\left(r^{\prime },\omega \right)·\mathit{d}\mathit{l}\end{array} \end{gathered}$$

(4)

where $\mathit{d}\mathit{l}$ is the wire path.

The first term in the last line of Equation (4) is the anomalous voltage:

$$\Delta V_j^a=\Delta V_{M_k{M}_l}^a\left(\omega \right)=\underset{M_k}{\overset{M_l}{{\displaystyle \int }}}\left(\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\stackrel{˭}{\mathit{G}}}_E\left(r^{\prime },r_i,\omega \right)\Delta \tilde{\sigma }\left(r_i,\omega \right)\mathit{E}\left(r_i,\omega \right)dV\right)·\mathit{d}\mathit{l}$$

(5)

and the second term is the background voltage:

$$\Delta V_j^b=\Delta V_{M_k{M}_l}^b\left(\omega \right)=\underset{M_k}{\overset{M_l}{{\displaystyle \int }}}{\mathit{E}}^{\mathit{b}}\left(r^{\prime },\omega \right)·\mathit{d}\mathit{l}$$

(6)

The line integrals in Equations (5) and (6) are evaluated using Gaussian quadrature. This creates a computational problem, as the number of different combinations of receivers is rather large, as each receiver can be referenced to every other receiver to create a data point. For the array shown in Figure 5, this leads to $\frac{N\left(N-1\right)}{2}=45$ unique voltage measurements where $N$ is the number of receiver electrode locations. This is repeated for each transmitter position for this array. One can imagine that this number grows very rapidly for large, distributed arrays.

This computational difficulty can be overcome by breaking the line integrals in Equations (5) and (6) into individual wire segments. Each segment of wire is defined by a straight line that is broken by either receiver electrodes or bends in the wire:

$$\Delta V^b{}_{M_k{M}_l}\left(\omega \right)={\displaystyle \sum }_{p=k}^{l-1}\underset{M_p}{\overset{M_{p+1}}{{\displaystyle \int }}}{\mathit{E}}^{\mathit{b}}\left(r^{\prime },\omega \right)·\mathit{d}\mathit{l}$$

(7)

where p is the dummy index of the receiver wire segment end points (electrodes or bends). In Equation (7), each wire segment is evaluated with Gaussian quadrature for each transmitter position.

$$\Delta V^b{}_{M_p{M}_{p+1}}\left(\omega \right)={\displaystyle \sum }_{q=1}^{N_G}{w}_q{\mathit{E}}^{\mathit{b}}\left(r_q,\omega \right)·\widehat{l}$$

(8)

where $w_q$ are the weights for the Gaussian quadrature, $r_q$ are the locations of the nodes, and $N_G$ is the order of the quadrature.

The number of Gaussian quadrature points is adjusted from three for simple dipole–dipole type configurations to 50 when the receiver wire and transmitter wire are parallel and close (maximum IC). For each unique combination of receiver positions, the voltage response is a summation of the appropriate wire segments. Every possible receiver pair can then be computed from a linear combination of the appropriate wire segments:

$$\Delta V^b{}_j\left(\omega \right)=\Delta V^b{}_{M_k{M}_l}\left(\omega \right)={\displaystyle \sum }_{p=k}^{l-l}\Delta V^b{}_{M_p{M}_{p+1}}\left(\omega \right)$$

(9)

In this way, excellent accuracy is maintained, yet the maximum number of segments is the number of receiver positions plus the number of bends in the wire. This is much less than the number of unique receiver pairs. This greatly reduces the computational load for the background field computation.

3.3. Numerical Modeling of Inductive Coupling in the Frequency and Time Domain

We begin with the simple case of a polarizable half-space for simplicity in developing some intuition about the IP and IC problems across a variety of source-receiver geometries and time and frequencies. The transmitter wire is modeled as a series of electric dipoles on the surface of the homogeneous earth and integrated with Gaussian quadrature. The receiver wires are modeled by integrating the electric field along the wire path from receiver electrode M to electrode N to synthesize the voltage difference which would be measured in the field:

$$\Delta V_{MN}={{\displaystyle \int }}_M^N\mathit{E}·\mathit{d}\mathit{l}$$

(10)

where electric field, E, is computed as discussed in Section 3.2 above.

We consider first a simple case of a half-space characterized by the Cole–Cole model parameters listed in

Table 1. Base Cole–Cole model parameters for half-space IC testing.

Parameter	Value
DC Resistivity (ρDC)	100 [Ohm·m]
Chargeability (η)	100 [mV/V]
Time Constant (τ)	0.1 [s]
Relaxation Factor (C)	0.5 [ ]

.

The current is 1 A. The time domain response is computed by using a cosine transform of the imaginary part of the frequency domain spectrum, which gives the step-off response. Note that, an infinite step-off is used for these calculations for simplicity; but using a 50% duty cycle repeating waveform of limited bandwidth allows us to draw the same conclusions, and the true waveform is used for the field data inversion. Also note that, the time range and frequency range shown here is greater than available commercial systems, but it is instructive to look at the expanded spectrum.

To isolate the effects of the inductive coupling and wire path, we have computed four different responses shown. One is the true response which includes all effects, as would be measured in the field. This includes the electromagnetic inductive response (IC), the true wire path, and the complex conductivity of the earth. The other three responses are a straight wire path from the M to N receiver electrode positions (ignoring the true wire path), setting $\frac{\partial B}{\partial t}=\nabla \times E=0$ (ignoring inductive coupling), and ignoring induced polarization (chargeability is 0; earth has a constant conductivity over frequency). Note that, in the case of the DC approximation, the wire path does not matter, as the fields are conservative, and the line integral in Equation (10) is path independent.

We show an example of the severe inductive coupling effects in Figure 6. This shows the frequency domain apparent resistivity and phase as functions of frequency, and time domain voltage and apparent chargeability. Time domain apparent chargeability is the voltage at some time after turnoff normalized by the on-time voltage. There are several features to point out in this figure. Panel (a) shows the geometry of the setup. Panel (b) shows that the DC approximation is excellent to around 1 Hz, at which IC begins to dominate. Panel (c) shows the frequency domain quadrature response. As one should expect, the phase reflects both the polarization and IC much stronger than the in-phase response. With this configuration, the IC dominates the phase for nearly the entire frequency range. Only at the very lowest frequencies (<3 × 10⁻⁴ Hz) does the IP affect the response. Above this frequency, IC effects completely dominate the spectrum. This frequency is well below the range of commercial systems. The same IC domination of the response is also seen in the space of frequency domain apparent resistivity (panel (d)) and phase (panel (e)). In the time domain (panels (f) and (g)) the IC part completely dominates from the earliest times to 10 ms. By 20 ms, the response is free of inductive coupling and only the IP decay remains.

Modeling shows that, as expected, longer offsets, receiver wire paths which traverse close to the transmitter dipoles, and more conductive media all increase the effects of the IC, both by increasing the amplitude and pushing the domination to later times or lower frequencies.

Figure 7 shows a similar configuration to Figure 6, but with the transmitter wire not running directly along the receiver wire. The responses share many characteristics of the response shown in Figure 6. In the frequency domain case (panels b–e), we can see that the induced polarization is important below 1 Hz, while above 3 Hz the IC dominates with little or no sensitivity to the IP. There is a small transition range where both IP and IC share contributions to the response spectrum. It is very interesting that, in the time domain response, the separation between the IP and no IP cases is complete and spans the entire time range of the field data collected (40 ms to 1.92 s). Beyond 20 ms, the true response, straight line approximation, and DC approximation all overlie each other. The only curve deviating from this behavior is the response which does not include IP effect. This means inductive coupling does not contaminate the IP data (for this configuration) after 20 ms. There is a small IC effect in the early times of the response, but it is relatively minor and disappears after a few milliseconds.

3.4. Discrimination of the IP Time Constant and Relaxation Coefficient When Considering the IC Effects

The goal of collecting the spectral IP data is in measuring the details about the polarization spectrum of the medium. Real rocks exhibit time constants that range for many orders of magnitude [17,18], while many time-domain surveys measure less than two decades of decay (e.g., 50 ms to 2 s). Additionally, inductive coupling can dominate the early times and higher frequencies, so this leads to the question if a reasonable range of time constants of the media can be determined from the real data.

Figure 8 shows the response of the half space with the parameters given in Table 1 with the same survey configuration as shown in Figure 6. Only the true wire path with full EM effects have been computed. Each curve shows a different time constant of the medium. These vary from 1 ms to 100 s, quite a wide and realistic range. Also shown is a case with no chargeability. Note that, for the frequency domain, there is no separation of the curves at frequencies above 20 mHz, which indicates that the time constant for this configuration would not be resolved in the frequency domain using conventional frequency ranges. In the time domain, there is separation of the curves for time constants from around 1 ms to 10 s. This occurs over the time window from 40 ms to 5 s. Outside of this range the time constant cannot be resolved. Longer time constants would require lower base frequencies and resolving shorter time constants would require accounting for the inductive coupling.

Figure 9 shows the response of the half space with the parameters given in Table 1 but with a 10 Ohm·m half-space. Again, the same survey configuration as shown in Figure 6 is used. For the frequency domain, there is no separation of the curves at frequencies above 2 mHz, which is a decade lower than when the earth has a resistivity of 100 Ohm·m. In the time domain, there is separation of the curves for time constants from around 10 ms to 10 s. This occurs over the time window from 300 ms to 5 s. This is much later than when the earth has a resistivity of 100 Ohm·m, and also is nearly a decade later in than the first time sample. In this case, inductive coupling would need to be accounted for to resolve short time constants and to make use of the earlier time channels.

Figure 10 examines the resolution of the relaxation parameter. The colored lines show the responses of half space with the properties shown in Table 1, but with C varying as shown in the legend. In the frequency domain, the effect of varying C is largely pushed outside of typical survey bandwidth, and the base frequencies would need to approach 1000 s to detect the variations. This is outside the patience of most field crews! The effect of increasing C from 0.1 to 0.5 in this figure is large enough to increase the magnitude of the response. Thus, over the limited time range available with most instruments, variations of C below 0.5 will be nearly indistinguishable from variations in chargeability. For large values of C (0.75 to 1.0), the inflection in the response becomes apparent (the dotted pink line), and this could be resolvable.

Thus, our study of the inductive coupling (IC) effect on the IP data shows that, the frequency-domain response has a much larger contribution from IC than the time-domain response. The time domain shows better sensitivity to IP and to the time constants. The wire path is less important in time domain modeling but can still have a large effect when transmitter receiver separations are large, and the medium is conductive (<100 Ohm·m).

In both frequency domain and time domain, there is about one decade of time or frequency response that contains information from both IC and IP. This means that by including wires paths and the full EM solution in forward modeling and inversion, we can extract IP information from the response one decade earlier than if the DC approximation is used. This is especially important when short time constants are expected.

Additionally, all the IP parameters are coupled, and this must be kept in mind when interpreting inversion results. The chargeability must be relatively large for the time constant and relaxation parameters to be resolved. The time constant is poorly resolved if the relaxation parameter is small (0.5 or less), and the relaxation parameter and chargeability can also be traded for an equivalent response when the relaxation parameter is small.

4. Inversion of the IP Data

4.1. Modeling of the Inductive Coupling

The modeling method to obtain the domain fields uses the integral equation method [30], and our formulation includes the inductive coupling naturally. The modeling is based on the quasi-stationary (not static or DC) Maxwell’s equations. The IP effect is completely represented by modeling the conductivity as a frequency dependent and complex term. The modeling is done in the frequency domain and the results are transformed to the time domain with a cosine transform. The sources are modeled as an integration of current dipoles along the length of the transmitter wires. The receivers are modeled by integrating the electric fields along the path of the receiver wire. This modeling therefore includes inductive and capacitive coupling from the transmitting wires to the ground, the ground to the active receiving wires, and the transmitter wires directly to the receiving wires. Because of the very high impedance, we assume there is no feedback from the active receiving wires to the earth and back to the receiving wires, or from the passive wires to the active receiving wires. All filters in the instruments are included by applying them to the frequency spectrum before transforming the result to the time domain.

4.2. Inversion

As we discussed above, in the case of studying the IP effect, the conductivity has to be treated as a complex and frequency dependent function, $\sigma =\sigma \left(\omega \right)$, which increases significantly the number of unknown parameters of the inversion [34,35]. We can reduce this number by approximating the conductivity distribution using a GEMTIP model (1), or its simplified version (2) described in Section 2 above. For this field data set, we invert for the GEMTIP parameters $\left({\sigma }_o,f,\tau ,C\right)$. These are all updated and included in the forward modeling during the inversion. At this point we are solving for the simplified GEMTIP model parameters, which is very similar to the Cole–Cole model. However, in future we plan on including elliptical parameters into the inversion, which provide a better representation of the complex resistivity of the rocks [16].

Our problem can now be written in the classic form of the operator equation:

$$d=A\left(\sigma f\left(\eta ,\tau ,C\right)-{\sigma }_b\right)=A_G\left(m\right)$$

(11)

where $A_G$ is GEMTIP forward modeling operator, and $m$ is a vector of the GEMTIP model parameters $\left[\sigma ,\eta ,\tau ,C\right]$.

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

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

(12)

where $W_d$ is the data weighting matrix, and $\alpha$ is a regularization parameter.

There are several choices for the stabilizer [28,29]. In this paper, for simplicity, we use the minimum norm stabilizer ($S_{MN}$), which is equal to the square L₂ norm of the difference between the current model $m$ and an appropriate a priori model $m_{apr}$:

$$S_{MN}\left(m\right)=\Vert W_m{\left(m-m_{apr}\right)\Vert }_{L_2}^2$$

(13)

where $W_m$ is the weighting matrix of the model parameters.

We used the re-weighted regularized conjugate gradient (RRCG) method to find the minimum of the parametric functional $P^{\alpha }$. The details of the inversion algorithm are described by [32,34]. This method is based on integral equation (IE) forward modeling method outlined above in Section 3.2, which requires the discretization of the volume containing the anomalous complex conductivity only. The background medium, surrounding the volume with the anomalous conductivity, is characterized by some known horizontally layered conductivity, ${\sigma }_b$. Another advantage of the IE method is that the sensitivities of the observed IP data to the changes of the anomalous conductivity can be expressed in close form using the electric and magnetic Green’s tensors defined for the background conductivity model ${\sigma }_b$.

Indeed, following [34] the sensitivity of the electric field at point j due to perturbations in the conductivity of the ith cell is given by the following formula:

$$\frac{\partial \mathit{E}\left(r_j,\omega \right)}{\partial {\sigma }_i}=\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\mathit{G}}_E\left(r^{\prime },r_i,\omega \right)\mathit{E}\left(r_i,\omega \right)dV.$$

(14)

The right hand-side terms contain the domain-to-receiver Green’s tensors (${\mathit{G}}_E$, 2nd rank) and the total electric field $\mathit{E}\left(r\right)$. This equation is the extended Born sensitivity. Do not confuse this with the Born sensitivity or Born approximation, which would use a 1D background field approximation instead of the total, rigorously calculated electric field used here.

The sensitivities for the voltage along this wire path are described, according to (6), by the following expression:

$$\frac{\partial \Delta V_{M_k{M}_l}\left(\omega \right)}{\partial {\sigma }_i}=\frac{\partial }{\partial {\sigma }_i}\underset{M_k}{\overset{M_l}{{\displaystyle \int }}}\mathit{E}\left(r^{\prime },\omega \right)·\mathit{d}\mathit{l}=\underset{M_k}{\overset{M_l}{{\displaystyle \int }}}\left(\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\stackrel{˭}{\mathit{G}}}_E\left(r^{\prime },r_i,\omega \right)\mathit{E}\left(r,\omega \right)dV\right)·\mathit{d}\mathit{l}.$$

(15)

The line integral in Equation (14) is evaluated using Gaussian quadrature.

$$\frac{\partial \Delta V_{M_k{M}_l}\left(\omega \right)}{\partial {\sigma }_i}={\displaystyle \sum }_{p=k}^{l-1}\underset{M_p}{\overset{M_{p+1}}{{\displaystyle \int }}}\left(\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\stackrel{˭}{\mathit{G}}}_E\left(r^{\prime },r_i,\omega \right)\mathit{E}\left(r_i,\omega \right)dV\right)·\mathit{d}\mathit{l}$$

(16)

where p is the dummy index of the receiver wire segment end points (electrodes or bends).

For the sensitivities, $\mathit{d}\mathit{l}$ remains constant over each integration because the receiver wires are assumed to be straight (or broken into straight segments), so the order of the tensor-vector products can be changed. This gives the following expression:

$$\frac{\partial \Delta V_{M_k{M}_l}\left(\omega \right)}{\partial {\sigma }_i}={\displaystyle \sum }_{p=k}^{l-1}\underset{M_p}{\overset{M_{p+1}}{{\displaystyle \int }}}\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\stackrel{˭}{\mathit{G}}}_E\left(r^{\prime },r_i,\omega \right)\mathit{d}\mathit{l}·\mathit{E}\left(r_i,\omega \right)dV.$$

(17)

Because the domain electric field is assumed to be constant over a cell, it can be removed from the domain integral. The domain electric field is also constant with respect to the line integral along the wire path and can be removed from integral as well:

$$\frac{\partial \Delta V_{M_k{M}_l}\left(\omega \right)}{\partial {\sigma }_i}={\displaystyle \sum }_{p=k}^{l-1}\left(\underset{M_p}{\overset{M_{p+1}}{{\displaystyle \int }}}\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\stackrel{˭}{\mathit{G}}}_E\left(r^{\prime },r_i,\omega \right)dV\mathit{d}\mathit{l}\right)·\mathit{E}\left(r_i,\omega \right).$$

(18)

The term in the parenthesis is the Green’s linear operator for the $i\mathrm{th}$ cell to the $p\mathrm{th}$ wire segment:

$${\mathit{G}}_{\mathit{i}}^{\mathit{p}}\left(\omega \right)=\underset{M_p}{\overset{M_{p+1}}{{\displaystyle \int }}}\underset{D_i}{\overset{}{{\displaystyle \iiint }}}{\stackrel{˭}{\mathit{G}}}_E\left(r\prime ,r_i,\omega \right)dV\mathit{d}\mathit{l}.$$

(19)

Note that, this is independent of both the total electric field in the domain, $\mathit{E}\left(r,\omega \right)$, and the source position. This means that this term only needs to be precomputed before the inversion and never updated, and that it only needs to be computed for each wire segment and frequency. The Green’s linear operator for each unique set of receiver positions can be computed by linear combinations of the wire segments, as is done for the background fields. The line integral in Equation (19) is computed through a Gaussian quadrature using 15 points per cell the wire passes over. Note also that this term is a first rank tensor (vector) because of the evaluation of the line integral in Equation (19).

The sensitivities are then composed as follows:

$$J_{ij}\left(\omega \right)=\frac{\partial \Delta V_{M_k{M}_l}\left(\omega \right)}{\partial {\sigma }_i}={\displaystyle \sum }_{p=k}^{l-1}{\mathit{G}}_i^p\left(\omega \right)·\mathit{E}\left(r_i,\omega \right)$$

(20)

and the anomalous voltage response can be found by combining Equation (19) with Equation (5):

$$\Delta V_j^a\left(\omega \right)={\displaystyle \sum }_{i=1}^{N_c}{J}_{ij}\left(\omega \right)\Delta {\tilde{\sigma }}_i$$

(21)

Equations (9), (20), and (21) provide the required quantities for modeling and inversion of the 3D distributed SIP survey data using voltages and exact wire paths in a computationally efficient manner.

The sensitivities for the GEMTIP parameters that we invert for $\left(\mathrm{e}.\mathrm{g}.,{\sigma }_o,f,\tau ,C\right),$ are found by the chain rule with Equations (1) or (2), as shown by [32,34].

The examples of synthetic IP model study and a case study of application of this method for the Copper Deposit in Mongolia were discussed in [32], where the reader can find additional details about the developed method of the IP data inversion.

5. Results of Interpretation of SIP Survey Data

We have applied the developed inversion code to these field SIP data. The inversion was run according to the description above with horizontal cell sizes of 50 m × 50 m. There were eight cells in the z direction varying from 26 to 223 m in thickness giving a total depth range of 0–600 m. There were nine transmitters and 23 receiver electrodes for total of 482 combinations and these each had 20 time channels from 40 ms to 1.92 s plus the on-time channel. The center line had the nine transmitters (see Figure 4) which were 200 m dipoles. The receivers were combinations of electrodes from all three receiver lines. The receiver lines were 350 m apart and approximately 3 km long.

The error model used for the observed data is given by the following equation:

$${\epsilon }_{\mathrm{i}}=\left|{\mathrm{d}}_{\mathrm{i}}\right|\ast 10\%+4E-6V/m$$

(22)

where $i$ indicates the $i\mathrm{th}$ data point or the error estimate of that point and $\left|{\mathrm{d}}_{\mathrm{i}}\right|$ is the absolute value of the ith data point. The data weights are the inverse values of these errors. The normalized ${\chi }^2$ at each iteration during the inversion is shown in Figure 11. The inversion converged from a normalized ${\chi }^2$ of around 3 to 1.1 at the 60th iteration, which was chosen to be optimal based on the noise level in the observed data, the flattening of the convergence curve, and examination of the resulting models. The optimal iteration is guided by theory but chosen by experience.

The field data collection was described in the previous sections. There is no standard plotting conversion for this type of distributed array data. We are plotting the physical fields to develop some intuition about the processes and to visually compare the observed and predicted data to analyze the data fit. Each image shows the voltage corresponding to one transmitter position. The plotting point is the average position of the two potential measurement points, as was shown in the synthetic example.

Figure 12 shows examples of the data fit for different transmitters and time channels. The observed and predicted data visually match very well. This is also demonstrated by the final ${\chi }^2$ of 1.1. A perfect fit would theoretically have a value of 1.0. Figure 13 presents examples of observed (solid lines) and predicted (dotted lines) voltage decays from the on-time measurement to the latest time at 1.8 s.

The inversion results under the transmitter line are shown in Figure 14. The inversion results of 3D distributed array data are shown in the top panel, the 3D inversion results from the 2D lines are given in the middle panel. The bottom panel shows the results of 3D from the Spectrum time domain airborne electromagnetic (TD AEM) data (inverted using the method of [36]). There were two boreholes drilled in the survey area. The position of one of these drill holes is shown in Figure 14 overlaid on the resistivity images. We do not have the coordinates of the other borehole.

All three inversion results show generally the same features. First, there is pervasive conductive overburden in this area. This is to be expected in this highly weathered environment. Most areas below this conductive overburden are resistive, but there are notable exceptions. The most interesting feature is a strong conductor near 700 mE. The depth of this conductor is different in the images produced by different datasets; however, the overall locations and sizes of the conductor are similar, which is remarkable considering that we compare the ground and airborne data. The top and bottom panels are the most similar in the top 150 m. A conductor is shown at depth in the middle panel, but note that the intersection of the borehole shows basalt with no pyrite at this location, which means this should probably be resistive. This is likely an off-line (3D) effect manifesting in the 2D results. All three methods image the conductive overburden. All three also show a strong resistor under the hill at 1400 mE. We expect the AEM to have a higher depth of investigation, and we have truncated the AEM results for this figure to 300 m for comparison. This greater depth of exploration explains the conductor at depth between 1200 and 1500 m in elevation which is only weakly imaged in both resistivity inversion results.

Figure 15 shows the chargeability results of the SIP inversions from the distributed array (top panel), and from the 2D dipole–dipole line data (bottom panel). The results have similarities, especially in the near surface from 1400 to 2000 mE, but they show differences in the strength of the chargeability and the 2D line shows the chargeable body much more broken. In the drilling location with pyrite confirmed by visual inspection of the core, the 2D dipole–dipole survey shows a small increase in chargeability while the distributed array clearly shows a significant increase of chargeability which also corresponds better to the drilling result. The large chargeability anomaly at 1500 mE has not yet been drilled due to logistical challenges for the location.

The fact that both the resistivity and the chargeability from the 3D IP inversion matches the drilling better than the 2D IP inversion demonstrates an importance of using the data collected by a distributed array for interpretation of the IP survey data and inverted in full 3D.

Regarding the time constant and relaxation parameter, these images contain very interesting features. The 3D inversion of the 2D line data (bottom panels) shows quite similar features and characteristics across the three parameters of chargeability, time constant, and relaxation parameter. This can be explained by the results in the synthetic section—C and chargeability are largely coupled and required low noise and accurate modeling to separate. Time constant is also difficult to separate from chargeability, especially when the relaxation parameter is 0.5 or less.

However, the 3D inversion of the distributed array data tells a different story. The four recovered parameters of the inversion all paint different pictures of the subsurface. We believe this is because they are all imaging different properties and provide independent information. The resolutions of the IP parameters are all linked and must be interpreted with this in mind. For example, neither the time constant nor the relaxation parameter has resolution when the chargeability is week. The time constant is also poorly resolved if the relaxation parameter is small (0.5 or less), and the relaxation parameter and chargeability are also coupled when the relaxation parameter is small.

Figure 14, Figure 15, Figure 16 and Figure 17 shows the general schematic of the geology encountered while drilling (borehole is shown in these figures overlying the geophysical results). The uppermost lithology is a weathered shale, which is moderately conductive and not chargeable. The basalt with quartz has lower resistivity (due to the pyrite) and is chargeable (due to the pyrite). Below this is basalt with little presence of pyrite. This account for its lack of chargeability. The basalt is moderately weathered which could explain an elevated conductivity as shown by the 2D lines, but more likely there is alteration or sulfides which are being imaged and were not intersected (off-line effects).

Figure 18 shows a picture of the quartz veining and associated sulfides minerals (pyrrhotite, pyrite, and chalcopyrite). The pyrite and chalcopyrite mineralization is responsible for the observed IP effect. The preliminary analysis of drilling results also shows quartz veins filled with the minerals, like pyrrhotite, pyrite, and chalcopyrite, with the presence of gold and silver. This points to a possible epithermal gold and silver vein type deposit. This is further encouraged by the proximity to the Mahd Ahd Dahab and Jabal Sayid deposits, aka the “Cradle of the Gold”. In the Mahd Ahd Dahab deposit, the polymetallic mineralization of chalcopyrite, sphalerite, galena, and gold and silver tellurides, occurs in the quartz veins as well. More drilling and rock sample analysis, however, will be required for this prospect area to make an estimate of the economical potential of the discovery.

6. Conclusions

In this paper, we have presented the results of the SIP survey conducted over one of the prospective targets located by the Glass Earth Pilot Project in Saudi Arabia. We have developed a novel 3D SIP inversion method and code which consider all inductive coupling, the wires paths, and includes IP effects through the use of the GEMTIP model. We have also demonstrated an importance of taking into account the inductive coupling (IC) effect in interpretation of the data acquired by the SIP survey with distributed array system with the synthetic modeling study.

Since the resistivities in this survey are generally greater than 100 Ohm·m and with the time constants on the order of 30 to 100 ms, the IP and IC responses did not overlap with the time windows and geometry of this survey. This means that in this specific case history, inductive coupling could have been small. This is not the case in general and this must be determined based on each specific survey’s geometry and expected physical properties. The important result of our study is that, it is impossible to know if the IC coupling affects the data until the survey has been completed and the data analyzed. Thus, it’s best to take the IC effect into account in analyzing the data. The inversion algorithm presented in this paper makes it possible to take into account the IC effect automatically, if it exists in the observed data.

The distributed systems have become widely used in industry, which makes this result of special practical significance. The 3D GEMTIP inversion with automatic accounting for the inductive coupling recovered a resistivity model which closely matched independently recovered resistivity models obtained previously by 3D inversion of the airborne EM data and matched more closely to drilling results than a 3D inversion algorithm of the 2D line data which ignoring wire paths and IC.

The Glass Earth Pilot Project has followed an exploration workflow from country wide selection through ground SIP survey, borehole drilling and discovery of blind mineralization by using the advanced 3D inversion methods detailed in this paper.

While the discovery is fortuitous, it would not have happened without the techniques of modern geophysics. There is no surface expression of this mineralization. More geologic analysis will need to be done in this area before we have a solid understanding of how exactly the time constant and relaxation parameter are related to geology, but the initial images are promising. The 300 m drill hole were drilled in an original area of nearly 8000 km². The drilling results identified quartz veins filled with the minerals, like pyrrhotite, pyrite, and chalcopyrite, with the presence of gold and silver type of mineralization in our survey area which is similar to an epithermal, low-sulfidation, polymetallic type deposit of Mahd Ahd Dahab Mine.