- Split View
-
Views
-
Cite
Cite
A Soueid Ahmed, A Revil, 3-D time-domain induced polarization tomography: a new approach based on a source current density formulation, Geophysical Journal International, Volume 213, Issue 1, April 2018, Pages 244–260, https://doi.org/10.1093/gji/ggx547
- Share Icon Share
Summary
Induced polarization (IP) of porous rocks can be associated with a secondary source current density, which is proportional to both the intrinsic chargeability and the primary (applied) current density. This gives the possibility of reformulating the time domain induced polarization (TDIP) problem as a time-dependent self-potential-type problem. This new approach implies a change of strategy regarding data acquisition and inversion, allowing major time savings for both. For inverting TDIP data, we first retrieve the electrical resistivity distribution. Then, we use this electrical resistivity distribution to reconstruct the primary current density during the injection/retrieval of the (primary) current between the current electrodes A and B. The time-lapse secondary source current density distribution is determined given the primary source current density and a distribution of chargeability (forward modelling step). The inverse problem is linear between the secondary voltages (measured at all the electrodes) and the computed secondary source current density. A kernel matrix relating the secondary observed voltages data to the source current density model is computed once (using the electrical conductivity distribution), and then used throughout the inversion process. This recovered source current density model is in turn used to estimate the time-dependent chargeability (normalized voltages) in each cell of the domain of interest. Assuming a Cole–Cole model for simplicity, we can reconstruct the 3-D distributions of the relaxation time τ and the Cole–Cole exponent c by fitting the intrinsic chargeability decay curve to a Cole–Cole relaxation model for each cell. Two simple cases are studied in details to explain this new approach. In the first case, we estimate the Cole–Cole parameters as well as the source current density field from a synthetic TDIP data set. Our approach is successfully able to reveal the presence of the anomaly and to invert its Cole–Cole parameters. In the second case, we perform a laboratory sandbox experiment in which we mix a volume of burning coal and sand. The algorithm is able to localize the burning coal both in terms of electrical conductivity and chargeability.
1 INTRODUCTION
In hydrogeophysics, the characterization of subsurface geological structures (geometry and petrophysical properties) is nowadays routinely performed by the means of geoelectrical methods such as Electrical resistivity tomography (ERT) and induced polarization (IP) techniques (e.g. Michot et al. 2003; Comte et al. 2010; Fiandaca et al.2012, 2013; Kemna et al. 2012; Binley et al. 2015). ERT is restricted to the mapping of the electrical resistivity field only. Electrical resistivity depends on many factors such as salinity, temperature, water content, and the cation exchange capacity (CEC) of the material (e.g. Waxman & Smits 1968; Shainberg et al.1980; Revil et al. 1998). While resistivity monitoring can be applied to the monitoring of the water content, for instance in agriculture (e.g. Michot et al. 2003), the interpretation of DC resistivity data alone is notoriously difficult. IP goes further by being able to map, in addition to the electrical resistivity, other physical parameters of interest such as the chargeability and a distribution of relaxation times (Kemna et al. 2012). In that sense, IP can be considered as a useful and fruitful extension of the conventional ERT.
IP effects were first discovered during the last century (e.g. Schlumberger 1920; Dakhnov 1962). It was observed that when injecting a primary current into the ground, and then suddenly shutting it off, porous soils and rocks are able to store reversibly electrical charges and produce a secondary voltage that is decaying over time. This secondary voltage can last for few seconds to few minutes depending on the duration of the impressed primary current and IP characteristics of the subsurface (Kemna et al. 2012). At the pore scale, this process is diffusive and the charge carriers go back to their equilibrium state driven by chemical potential gradients once they have accumulated at some polarization length scales (e.g. at the edges of grains or pores).
Historically, the IP method was mainly used in mineral exploration for the detection of ore deposits because the chargeability of these targets is generally very strong (e.g. Bleil 1953; Van Voorhis et al. 1973; Zonge & Wynn 1975; Pelton et al. 1978; Telford et al. 1990). Later on, thanks to the technological progresses made in data acquisition, sensitivity of the instruments, and computers performances (e.g. Zimmermann et al. 2008), IP has become a very important method to investigate a broad spectrum of environmental applications. One can cite, for instance, the study of contaminants plumes (e.g. Olhoeft 1984, 1985, 1986; Vanhala et al. 1992, Vanhala 1997; Slater & Lesmes 2002; Kemna et al. 2004; Wainwright et al. 2014), the characterization of permeability and pore size distribution (e.g. Sturrock et al.1999; Revil et al. 2015c; Joseph et al.2016; Osterman et al. 2016), and recently coal seam fires detection and localization (e.g. Shao et al. 2017) just to cite few examples among a very rich literature.
The conventional time-domain induced polarization (TDIP) is restricted to the evaluation of the DC electrical conductivity σ0 and the chargeability M. However, the IP properties of soils and rocks can be represented by a distribution of relaxation times as well. This distribution can be sometimes simplified by a mean and a standard deviation. For instance, in a classical representation model known as the Cole–Cole model (Cole & Cole 1941), the distribution of relaxation times is described by two parameters namely the relaxation time τ, which describes a mean relaxation time, and the Cole–Cole exponent c, which describes the broadness of the distribution. The relaxation timeτrefers to the main time taken by a material that has been submitted to an electrical field or an electrical current, to go back to its equilibrium state. In a Debye model the relaxation time is the time required to see the secondary voltage falling down by a factor exp(−1) from its nominal value. Thanks to physical models such as the dynamic Stern layer model (Rosen & Saville 1991; Rosen et al.1993; Revil 2012, 2013), these Cole–Cole parameters can be interpreted in terms of textural and electrochemical properties of the material under consideration. The dynamic Stern layer model has proven indeed to be an efficient model to explain various empirical trends observed in the literature for rocks in absence of metallic particles (see for instance the works by Weller et al.2011, 2013, 2015a,b, who developed a series of empirical relationships all explainable by the dynamic Stern layer model). In the absence of metallic particles, this includes a mean pore (grain) size, a pore (grain) size distribution and the CEC of the material, which properties can be independently measured in the laboratory (experimental checks are for instance provided by Revil et al. 2014; Niu & Revil 2016). Induced polarization can also bring information regarding the presence of semi-conductors, metals, and semi-metals (Pelton et al. 1978; Revil et al. 2015a,b; Mao & Revil 2016; Mao et al. 2016; Revil et al. 2017a,b). This is due to the very strong IP response associated with the presence of metallic particles embedded in a porous material, which can be also affected by redox processes (Wong 1979) and the polarization of the pore water around the metallic particles (Misra et al. 2016a,b).
A number of published works have been conducted to image the Cole–Cole parameters in the subsurface. For instance, Loke et al. (2006) used a 2-D least square inversion to recover the Cole–Cole parameters distributions in a laboratory sandbox experiment. Ghorbani et al. (2007) inverted the Cole–Cole parameters using a 1-D Bayesian inference approach and they applied their methodology on synthetic homogenous half spaces. Yuval & Oldenburg (1997) estimated, in 2-D, the Cole–Cole parameters from TDIP data using a very fast simulated annealing approach. They successfully recovered these parameters on synthetic and real field data sets. Fiandaca et al. (2012) developed a forward and inverse code, which takes into account the modelling of the transmitter waveform and the receiver transfer function. Such methodology allowed for improving the resolution of the estimated Cole–Cole parameters. Recently, Nivorlis et al. (2017) proposed a 3-D computation scheme to retrieve in 3-D the Cole–Cole parameters. The last step of their work is accomplished using a particle swarm optimization algorithm. All these works follow the same path in describing induced polarization in terms of a time dependent electrical conductivity problem. In this approach, conductivity changes from an instantaneous conductivity (all the charge carriers are mobiles) to a DC conductivity for which some of the charge carriers are blocked at some polarization lengths scales and do not participate anymore to the conduction process. This approach finds its roots in the seminal work of Siegel (1959), who proposed to model the IP effects as a perturbation of the electrical conductivity field by the chargeability. Following this approach, an apparent chargeability is obtained by solving the Poisson's equation twice: once with the DC electrical conductivity σ∞ as input (σ∞ denotes the instantaneous conductivity, induction effects being neglected) and the other by taking σ0 = σ∞(1 − M) as input (i.e. using the DC conductivity distribution). This method has the merit of being straightforward and uses the same forward operator corresponding to the Poisson equation to solve the conductivity and chargeability problems. This formulation has been widely taken up and used by a majority of geophysicists.
Our approach follows a quite different path, which can also be traced to the seminal work of Siegel (1959). Indeed, Siegel demonstrated that the (primary) current injection creates a secondary current density JS(t) in the conductive ground. This secondary current is related to the primary current Jp through the chargeability evolution once the primary current has been shut off. This point, first raised by Siegel (1959) to our knowledge, has not been used by IP practitioners (despite the advantage that comes with it as discussed below). Since the secondary source current density is formally similar to a diffusion source current density in self-potential studies, time domain induced polarization can be described as a time-dependent equivalent self-potential-type problem (see Mao & Revil 2016, for a preliminary step in this direction). We can be more explicit here. The secondary source current density is driven in induced polarization by chemical potential gradients exactly like diffusion potentials in self-potential studies (e.g. Ikard et al. 2012). The only formal difference is that in induced polarization, the chemical potential gradient of the charge carriers have been ‘actively’ set up by the injection of the primary current (through cross-coupling effects, see for instance Revil 2017a,b), while in classical self-potential studies, the chemical potential gradients can come from the injection of a salt tracer in the environment (e.g. Jardani et al. 2013).
The advantage of formulating the IP problem as an equivalent self-potential problem is that we avoid solving a nonlinear inverse problem because retrieving the current density JS(t) from the recorded electrical field is a linear problem. Examples of such linear problem in self-potential tomography can be found in Mahardika et al. (2012) and Haas et al. (2013) and in electroencephalography for instance by Trujillo-Barreto et al. (2004). A linear inverse problem does not require the use of an iterative process. This means that notable computational time savings can be made for the tomography of the intrinsic chargeability field. New instrumentations that are massively multichannels such as the IRIS Full waver instrument can operate the way we advocate: all the stations measure simultaneously the electrical field (i.e. the gradient of the electrical potential distribution along the curvilinear coordinates of the ground surface) for each injection bipole [A, B].
In this work, the IP data collection is performed in a self-potential ‘fashion’, this means that a limited number of primary current injections is performed and the secondary voltage measurements are recorded at all remaining electrodes like in a self-potential survey can save a lot of time with modern multi-channel technologies. Considering the secondary voltages as pseudo self-potential data is also correct from a physical point of view since these secondary voltages are driven by chemical potential gradients. They are therefore identical in nature to diffusion potentials as mentioned above.
The goal of this paper is to develop the novel approach mentioned above by formulating the TDIP forward and inverse problems as an equivalent self-potential problem. We present a 3-D framework for recovering the Cole–Cole parameters spatial distributions from TDIP data. To the best of our knowledge, this is the first attempt to recover the Cole–Cole parameters in 3-D using such an approach. The proposed algorithm is validated on two cases. (i) A synthetic model where the Cole–Cole parameters true distributions are known and will be compared to the estimated ones. (ii). A laboratory experiment is performed with some coal burning in a sandbox. Our goal is to develop a proof-of-concept of the method before to explore complex geometries in future contributions and to be didactic in describing the step-by-step procedure in getting the end-results.
2 FORWARD MODELLING
2.1 Electrical conductivity and chargeability
Once eq. (1) is solved, one can compute the apparent resistivity associated with the instantaneous conductivity through ρa = KψMN/I, where K is the geometric factor (depending on the position of A, B, M, and N and the topography), ψMN is the potential difference between two measuring electrodes M and N (Fig. 1a). Furthermore, eq. (1) can be seen as a nonlinear mapping operator F(.) associating the electrical potential |${\psi _{{\sigma _\infty }}}$| (Fig. 1) to the electrical conductivity σ∞, that is, |${\psi _{{\sigma _\infty }}} = F( {{\sigma _\infty }} )$|. This potential |${\psi _{{\sigma _\infty }}}$| is instantaneously recorded when the current injection is turned on. Similarly, a potential |${\psi _{{\sigma _0}}} = F( {{\sigma _0}} )$| can be registered when the primary current has been applied long enough (Fig. 1). The concepts of instantaneous and DC conductivities are explained physically in Fig. 1(b) in the context of the dynamic Stern layer model.
Therefore, in the classical approach, modelling the time domain IP requires solving the electrical conductivity equation (1) twice with two distinct electrical resistivity distributions (as proposed by Siegel 1959).
Note that this approximation is valid only under the condition that ti, ti + 1 ≪ τ (see Florsch et al.2011). This equation can be used to determine for each quadripole ABMN (A and B being the current electrodes and M and N the potential electrodes), an apparent chargeability Ma. Note that our notations are pretty standard. Some authors used sometimes the letters m or η to denote the chargeability.
2.2 Forward modelling in the charging phase
Eq. (11) means that the dipole moment associated with the polarization of a grain (see Fig. 1b) is antiparallel to the applied current density, explaining the sign ‘−’ in this equation. Further details on eq. (11) can be found in the work of Seigel (1949, 1959) and will not be repeated here. Note that we are assuming a Cole–Cole complex resistivity model below which, in turn, can be easily related to a Cole–Cole complex conductivity model using a relationship between the time constants (see Florsch et al.2012; Tarasov & Titov 2013). Eq. (11) is the key equation of this paper. The source (secondary) current density (index s) is intrinsically associated with cross-coupling phenomena associated with the existence of ionic chemical potential gradients at the scale of the grains or the pores (see Fig. 1b).
In the next section, we study the material response in the discharging phase.
2.3 Forward modelling in the discharging phase
3 COLE–COLE PARAMETERS ESTIMATION
The model of conductivity logarithm si is updated until the convergence criteria are met. Usually few iterations (in our experience ∼5 to 7) are required for reaching convergence.
The computation of the kernel is done as follows. For each cell (numbered from 1 to m), we assign in Comsol Multiphysics (using the finite elements method) an elementary dipole in the three orthonormal directions (x, y, z) (so three elementary dipoles in total) and we compute the resulting distributions of the potential at each of the nϕ voltage electrodes recording the secondary voltages. In each case, we remove the potential at the position of the selected reference electrode (i.e. the reference station used to reference the secondary voltage for the nϕ − 1 scanning electrodes). Note that the total number of electrodes will be nϕ + 2 accounting for the two electrodes used to inject the primary current and generally not used to measure the secondary voltages. As explained in details in Jardani et al. (2008), the computed kernel should respect the potential at the selected reference electrode (i.e. equal to zero at all times). The kernel is composed of three matrices G = [Gx, Gy, Gz] each of these matrices Gi(i = x, y, z) is a nϕ × m matrix so G corresponds to a N × 3 M matrix. The sources will be described by the current dipole moment vector p = ID where D denotes the displacement vector pointing in the direction of the flow of the current (in the direction of the electrical field) and p the current dipole moment vector (this is equivalent to the current density vector divided by the volume of the cell). The current dipole moment is therefore expressed in A m. An important mark is that this step is the most tedious but is performed only once as soon as the conductivity distribution has been determined.
Once the distribution of JS(t) has been recovered, we divide it by the distribution of the primary current Jp to obtain the intrinsic chargeability M(t). Along these lines, it is important to point out that, dividing JS(t) by Jp may lead to severe artefacts if the values of Jp are too small. This may be the case in the region that are far from the sources and where the sensitivity is too low. In practice, we set a critical value of Jp and, we neglect all the values smaller than this threshold. Once the intrinsic chargeability decay curve is obtained on each cell, we can perform the estimation of the remaining Cole–Cole parameters, that is, M0, τ and c. We opt for the strategy of using a least-square curve fitting of the intrinsic chargeability decay (which follows a Cole–Cole relaxation model) on each cell, to obtain a value of M0, τ, c.
Fitting the chargeability curve can lead to values of M0, τ and c that are not physically meaningful but still, perfectly match the intrinsic chargeability curves decay. To overcome this issue, it is pertinent to impose constraints on M0, τ and c. We have 0 < c < 1, 0 < M0 < 1 and τ > 0, but in practice, generally c ∈ [0, 1] in the presence of metallic particles and c ∈ [0, 0.5[ in the absence of those particles. In our case, we imposed on each cell j the constraints as lower and higher bounds: M0, j ∈ [0, 1], cj ∈ [0.5, 1] and τj ∈ [0, + ∞]. The choice of such range for c is motivated by the fact that, similarly to metallic particles, the burning coal exhibits high polarization effects as we will see later on in the current work. We use a least square technique for this third step but we note that many algorithms have been designed in the literature to solve this problem (e.g. Ghorbani et al.2007; Chen et al. 2008; Keery et al. 2012; Bérubé et al. 2017).
In summary, our approach is based on solving three inverse problems in cascade (Fig. 2). We first compute the conductivity field from which we compute the kernel for the secondary voltages. Then, we solve a set of inverse problems for each bipole current injection to determine the time-dependent chargeability of each cell. For each cell, we use a third inversion algorithm to retrieve the Cole–Cole parameters.
4 APPLICATION TO A SYNTHETIC CASE
We now present two cases on which we apply our approach. The first case corresponds to a synthetic test and we estimate the Cole–Cole parameters from the TDIP data, which are simulated using the forward modelling presented in Section 2. This synthetic test is convenient for benchmarking our method as the true fields are known and thus the comparison between these fields and the estimated ones can be readily done.
In this synthetic test, we simulate a sandbox with insulating boundaries. The central portion of the simulated tank is occupied by a polarizable body. We use a network of 24 electrodes placed around this body (Fig. 3) and the box is discretized with 1000 elements (10 × 10 × 10 elements). This electrode configuration has been used such as the primary current flow across the target and therefore allows for adequately discriminating its properties. The experiment data configuration is summarized in Table 1. The true Cole–Cole M, σ, τ and c distributions are provided in Fig. 4(a). These fields have been generated to create a contrast between the Cole–Cole parameters of the anomaly and those of the background, which is also polarizable but with a weaker polarization. A total of two current injections (described in Fig. 3) have been performed for the IP data acquisition. A current of 1 mA is injected for 2 s, then the secondary voltage decay is recorded on 8 time windows of 1 s each (see Table 2).
Injected current amplitude . | Number of injections . | Number of electrodes . | Delay time . | Injection duration . | Measurement times . |
---|---|---|---|---|---|
1 mA | 2 | 24 | 0.1 s | 2 s | 1 s, 2 s, 3 s, 4 s, 5 s, 6 s, 7 s, 8s |
Injected current amplitude . | Number of injections . | Number of electrodes . | Delay time . | Injection duration . | Measurement times . |
---|---|---|---|---|---|
1 mA | 2 | 24 | 0.1 s | 2 s | 1 s, 2 s, 3 s, 4 s, 5 s, 6 s, 7 s, 8s |
Injected current amplitude . | Number of injections . | Number of electrodes . | Delay time . | Injection duration . | Measurement times . |
---|---|---|---|---|---|
1 mA | 2 | 24 | 0.1 s | 2 s | 1 s, 2 s, 3 s, 4 s, 5 s, 6 s, 7 s, 8s |
Injected current amplitude . | Number of injections . | Number of electrodes . | Delay time . | Injection duration . | Measurement times . |
---|---|---|---|---|---|
1 mA | 2 | 24 | 0.1 s | 2 s | 1 s, 2 s, 3 s, 4 s, 5 s, 6 s, 7 s, 8s |
Injected current amplitude . | Number of injections . | Number of electrodes . | Delay time . | Injection duration . | Measurement times . |
---|---|---|---|---|---|
5 mA | 2 | 59 | 0.03s | 1s | 0.05 s, 0.09 s, 0.17 s, 0.33 s, 0.65 s, 1.29 s, 2.57 s, 5.13 s. |
Injected current amplitude . | Number of injections . | Number of electrodes . | Delay time . | Injection duration . | Measurement times . |
---|---|---|---|---|---|
5 mA | 2 | 59 | 0.03s | 1s | 0.05 s, 0.09 s, 0.17 s, 0.33 s, 0.65 s, 1.29 s, 2.57 s, 5.13 s. |
Injected current amplitude . | Number of injections . | Number of electrodes . | Delay time . | Injection duration . | Measurement times . |
---|---|---|---|---|---|
5 mA | 2 | 59 | 0.03s | 1s | 0.05 s, 0.09 s, 0.17 s, 0.33 s, 0.65 s, 1.29 s, 2.57 s, 5.13 s. |
Injected current amplitude . | Number of injections . | Number of electrodes . | Delay time . | Injection duration . | Measurement times . |
---|---|---|---|---|---|
5 mA | 2 | 59 | 0.03s | 1s | 0.05 s, 0.09 s, 0.17 s, 0.33 s, 0.65 s, 1.29 s, 2.57 s, 5.13 s. |
As illustrated on Fig. 2, we start by the key step of retrieving the spatial distribution of the conductivity field. This is a classical ERT problem in which the measured resistances are used as the observed data. Fig. 4(c) shows that the retrieved conductivity field reproduces very clearly the true conductivity field (both the geometry and the magnitude of the conductivity). Next, we use the conductivity field to compute the secondary (source) current density field Js. This is done via a linear inversion process as described in Section 3 and Fig. 2. The retrieved Js field is shown in Fig. 5(a), an area of high current density coincides with the localization of the anomalous body. The intrinsic time-dependent chargeability field is then obtained using eq. (41). While doing so, we pay attention to the behaviour of the primary current field Jp, because it strongly decays as we move away from the current electrodes, and therefore using eq. (41) out of the limits of the high sensitivity pattern depicted in the Jp field (see Fig. 5b) can lead to some artefacts. To cope with this issue, we simply discard the cells of Jp which are too far from the sensitivity pattern and which do not bring any relevant information about the localization of the target. The intrinsic chargeability M(t) is represented on Fig. 6. We observe that the anomaly area has a higher chargeability (as expected) and it keeps weakening as the polarization phenomenon evolves with time, to completely vanish at the last measurement times.
The final step of the inversion scheme is to jointly estimate the remaining Cole–Cole parameters, that is, M0, τ and c. We perform this step by doing a least-square constrained curve fitting of the intrinsic chargeability decay. The constraints that we use are non-restrictive and simply reflect the physical behaviour of these parameters. We impose that M0 ∈ [0, 1], τ ∈ [0, + ∞] and c ∈ [0, 1]. We represent in Figs 4(d)–(f) the results of such estimation, as one can see, the anomaly is remarkably well-recovered and the magnitudes match those of the true distributions.
This synthetic study case was indeed of high importance because it gave us the possibility to check the validity of our approach, for the forward problem as well as the inverse one. In the light the results of the synthetic experiment, it is quite logical to expect that our approach can be used for efficiently characterizing localized target that have high polarization effects. This point will be explored in more details on a sandbox experiment data set in the next section.
5 APPLICATION TO A SANDBOX EXPERIMENT
5.1 Experiment setup
The experiment consists in mixing two materials, a clean sand and some burning coal. The latter is highly polarizable (Shao et al. 2017), while the sand is weakly polarizable because characterized by a low CEC (Mao et al. 2016). Such high polarizability contrast is suitable for validating our procedure. We use a 46 cm × 29 cm × 27 cm plastic tank (with insulating boundary conditions on all the faces). We place a network of 52 stainless steel electrodes (so generally nϕ = 49, 52 electrodes less the reference electrode and two current electrodes). These electrodes are placed horizontally and vertically so they surround the central area of the sandbox (Fig. 7). Electrode #1 is used as a reference electrode for all the secondary voltage measurements. The box is discretized with m = 1000 elements (10 × 10 × 10 elements).
In the central area of the sandbox, a cylindrical volume placed 6 cm above the bottom boundary of the tank and extending over 12 cm of diameter and 15 cm of height is filled with coal chunks that are ignited before the induced polarization measurements are collected. We keep blowing air through the coal (using a hair dryer) until its combustion has well progressed. From now on, we will refer to this burning coal area as the target. In the aftermath, the target is delicately covered with moistened sand. We use silica sand that is composed of 95 per cent SiO2, 4 per cent KSi3AlO8 and less than 1 per cent NaAlSi3O8, is characterized by a mean grain diameter and standard deviation 130 ± 20 μm (porosity 0.34 and a formation factor of 3.6). The tap water used for humidifying the sand has an electrical resistivity of 18.9 Ohm.m at 25°C. This water gives the saturated sand a resistivity of around 68 Ohm.m. The reason behind humidifying the sand is that we need to make it more conductive in order to foster the electrical conduction through the tank.
We used an ABEM resistivity meter to acquire the IP data (apparent resistivity, apparent chargeability, secondary voltage). These IP data were acquired like in a self-potential survey, that is, we inject current between two electrodes A and B and we measure the potential between all the remaining electrodes with respect to the reference electrode (electrode #1). By doing so, we notably reduce the acquisition time as we only need a limited number of current injections. In our experiment, we only performed two dipole current tests using the current bipoles [A, B] = [5, 14] and [31, 40]. Thus, a number of 147 measurements were acquired in barely 5 min. This rapid acquisition protocol is particularly suitable for this kind of coal experiment, because the coal combustion is fast and measurements must be recorded before the extinguishment of the burning coal. For the 2 primary current injections, we inject a current of 5 mAduring 1 s, then, we shut off the current and we start recording the secondary voltage decay after a delay time of 0.13 s. This delay time is used to avoid spurious electromagnetic effects, especially capacitive and inductive coupling effects. Then the secondary voltage decay is measured using ten time windows of different lengths (see Table 2).
As for self-potential data (e.g. Jardani et al. 2008; Haas et al. 2013), the secondary potential voltages need to be referenced in time in order to avoid static electrical potential differences between the scanning electrodes and the reference electrode. At the last time window, we assume that the voltage decay has fully relaxed to reach zero, and we use it as temporal reference. In other words, we subtract the secondary voltage of the last time window (that is the voltages at 5.13 s) from all the secondary voltages measured during the other time windows. Figs 8(a) and 8(b) portray the distribution of the secondary voltages generated for the following current injections: [A, B] = [5, 14] and [A, B] = [31, 40], respectively. The burning coal shows a high polarization level compared to the sand in agreement with the recent study by Shao et al. (2017).
5.2 Inversion
We discretize the simulation domain into 10 × 10 × 10 rectangular cells and we apply no flux boundary condition n · ∇ψ = 0 (n is the outward unit vector at the boundary of the sandbox and ψ denotes the electrical potential) at all the borders of the sandbox. We estimate for each cell a value of the electrical conductivity, the source current density, the time dependent intrinsic chargeability, the relaxation time, and the Cole–Cole frequency exponent. At the end, the assembly of the values on each cell will give the fields of the corresponding parameters on the whole simulation domain.
For the electrical tomography survey (ERT), we inject the current on the electrodes [A, B] = [5, 14], and we measure the resistance on all the remaining electrodes, then we switch the current electrodes to the dipole [A, B] = [31, 40] and we perform the resistance measurements on the remaining electrodes. As described in Section 2, we use this observed resistance data to recover the electrical conductivity spatial heterogeneities. The inversion results are illustrated on Fig. 9. It shows an anomaly of high electrical conductivity (∼0.12 S m− 1) in agreement with the area where the coal is burning. In addition, the observed and computed resistances agree with each other (Fig. 10). This result indicates that when we are dealing with a localized target, such as in our experiment, only few pairs of current electrodes are sufficient to retrieve a satisfactory tomogram provided that the current flows through the target. In addition, examining the primary current distributions for the two bipole injections (see Fig. 11) shows a high sensitivity region going from one current injection toward the other one and passing through the target area.
Since JS(t) is now known and Jp is already given by Ohm's law as Jp = −σ∇ψ (E = −∇ψ denotes the applied electrical field) no additional inversion steps are required to retrieve the intrinsic chargeability. That said, special care must be taken to adequately perform the cell-by-cell division in eq. (49). In fact, as already illustrated in Fig. 11, as we get far from the current electrodes A and B, we lose the sensitivity and some cells away from the sensitivity pattern may have small values of Jp resulting though in some artefacts in the tomogram of M(t). The solution to this problem is to use a critical value of Jp below which, the corresponding primary current cell is not taken into account in the inversion for a given bipole [A, B].
We present in Figs 14(a)–(h) the intrinsic chargeability distributions at all the measurements times. We can observe that the high chargeability anomaly, which is likely associated with the burning coal has the highest magnitude just after the current was shut down when the polarization of the burning coal is the greatest, and keeps decaying to zero when the potential has fully relaxed during the last time window at 5.13 s.
At this stage, three parameters are still missing in our study. There are the relaxation time, the Cole–Cole exponent, and the chargeability at time zero, M0, which could not be estimated using the previous approach because we do not record the voltage measurement at time zero since we consider a delay time (i.e. the dead time contaminated by electromagnetic effects). As we have now computed the intrinsic chargeability decay curve on each cell, these three parameters can be estimated altogether by simply fitting the time domain Cole–Cole model in eq. (35). However, we impose some constraints as lower and higher bounds for each of these parameters (see Section 2). We used a least-square curve fitting in which M0, τ and c values on each cell were initially set to 0, 1 s and 0.5, respectively. The results of such estimation are represented on Fig. 15. We are only interested by the localization of the highly polarizable target, therefore, we only represent the cells that have a critical chargeability value which (>0.05). We note that the M0 field has higher values than the chargeabilities of all other times which, indicating that the coal is the most chargeable just after the current is shut off. The τ- and c-distributions are represented in the region characterized by the high values of the chargeability corresponding to the target. The values of c are close to 1 in the target area which suggests that a Debye model could be used to fit the intrinsic chargeability decay curve for coal. This result is in line with those obtained by Shao et al. (2017).
Plotting the intrinsic chargeability curves and the fitted model (see Fig. 16) shows that there is an excellent match between both, indicating the reliability of our inversion results. The relaxation time τ is comprised between 0.4 s and 1.4 s. Further investigations should be conducted to better understand the values of the relaxation time in terms of physical process and the nature of the charge carriers.
6 CONCLUSION
We have developed a novel approach for solving the forward and inverse time-domain induced polarization problem. Our methodology relies on the fact that the IP phenomenon is similar to a time dependent (highly transient) self-potential phenomenon including at its physics-based level. This brings out the rationale for reformulating the time-domain induced polarization problem as an equivalent self-potential problem. This methodology has notable advantages over the conventional way of considering the induced polarization effects: (i) We can drastically reduce the acquisition time by using a limited number of current injections and using all the other electrodes as potential electrodes sampling the secondary potentials with respect to a reference electrode. (ii) The computational time required for modelling the induced polarization effects is highly shortened. (iii) Using the present approach, we obtain a time dependent intrinsic chargeability distribution which allows for monitoring in real time, the evolution of polarizable targets. We successfully validated our approach by recovering the 3-D spatial heterogeneities of the Cole–Cole parameters in a numerical experiment and a laboratory sandbox experiment. The next steps will be to apply our methodology to complex geometries such as found in field conditions, to parallelize the procedure, and to test and compare various strategies to retrieve the Cole–Cole parameters at each cell.
Acknowledgements
This contribution is supported by the University of Melbourne through a project funded by the Commonwealth of Australia (contract CR-2016-UNIV.MELBOURNE-147672-UMR5275). We thank the Editor, Nicolas Florsch and an anonymous referee for their constructive and fruitful comments that we have used to make this manuscript more didactic.
REFERENCES