TY - JOUR
AU - Revil, A
AB - 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. Electrical properties, Hydrogeophysics, Electrical resistivity tomography (ERT), Inverse theory, Numerical modelling, Tomography 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 We consider below time scale and length scales such as the induction effect can be neglected. The primary electrical potential field generated just after a current injection in a 3-D heterogeneous isotropic medium can be described by the Poisson equation as   \begin{equation}\nabla \cdot \left( {{\sigma _\infty }\nabla \psi } \right) + I{\delta _3}({\boldsymbol{r}} - {{\boldsymbol{r}}_0}) = 0,\end{equation} (1)where σ∞(in  S m-1) denotes the 3-D high frequency electrical conductivity field (i.e. the instantaneous conductivity of a medium submitted to an electrical field, electromagnetic induction effects neglected, see Fig. 1), ψ (in V) is the electrical potential field generated by the injection of the current I (in A), δ3(r − r0) = δ(x − x0) δ(y − y0) δ(z − z0) is the Dirac distribution, x, y, z represent the space locations and x0, y0, z0 are the spatial coordinates of the current sources locations. Eq. (1) is subject to the following boundary conditions   \begin{equation}\psi = \alpha \quad {\rm{on}}\,{\Gamma _1},\end{equation} (2)  \begin{equation}{\sigma _\infty }\nabla \psi \cdot \,{\bf n} = \beta \quad {\rm{on}}\,{\Gamma _2},\end{equation} (3)with Γ1∪Γ2 = ∂Ω where ∂Ω denotes the simulation domain boundaries, n is the outward unit vector perpendicular to Γ2. In case of α = 0 and β = 0, we refer to these boundary conditions as homogenous Dirichlet and homogenous Neumann boundary conditions, respectively. If the subsurface is going to infinity, we can take the potential going to zero at the external boundary of the domain (this is easily done with a finite elements solver using coarse meshing outside the area of interest and performing some benchmark testing). Figure 1. View largeDownload slide Polarization of a porous material. (a) Recorded voltage at the two voltage electrodes M and N as a function of time. The primary current is applied for the period Ton (from –T to 0). The secondary voltage is measured after the primary current is shut down (during Toff). (b) Polarization of the grains. The instantaneous conductivity σ∞ is defined right after the application of the primary current. All the charge carriers are mobile. After a long time, the some of the charge carriers are blocked. This defined the direct current conductivity σ0. Figure 1. View largeDownload slide Polarization of a porous material. (a) Recorded voltage at the two voltage electrodes M and N as a function of time. The primary current is applied for the period Ton (from –T to 0). The secondary voltage is measured after the primary current is shut down (during Toff). (b) Polarization of the grains. The instantaneous conductivity σ∞ is defined right after the application of the primary current. All the charge carriers are mobile. After a long time, the some of the charge carriers are blocked. This defined the direct current conductivity σ0. 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. The chargeability distribution is defined as   \begin{equation}M = 1 - \frac{{{\sigma _0}}}{{{\sigma _\infty }}} = \frac{{{\sigma _\infty } - {\sigma _0}}}{{{\sigma _\infty }}},\end{equation} (4)which is equivalent to   \begin{equation}{\sigma _0} = {\sigma _\infty }\left( {1 - M} \right).\end{equation} (5) This equation has a clear physical meaning to the light of Fig. 1(b): in steady-state conditions, some of the charge carriers are blocked and the DC conductivity is reduced by a factor (1 − M) with respect to the instantaneous conductivity for which all the charge carriers are mobile. This reduction is also thoroughly discussed by Seigel (1949, 1959). The DC conductivity is obtained with   \begin{equation}{\psi _{{\sigma _0}}} = F\left( {{\sigma _\infty }\left( {1 - M} \right)} \right).\end{equation} (6) Similarly, an apparent chargeability can be computed as   \begin{equation}{M_a} = \frac{{{\psi _{{\sigma _\infty }}} - {\psi _{{\sigma _0}}}}}{{{\psi _{{\sigma _\infty }}}}} = \frac{{F\left( {{\sigma _\infty }} \right) - F\left( {{\sigma _\infty }\left( {1 - M} \right)} \right)}}{{F\left( {{\sigma _\infty }} \right)}}.\end{equation} (7) 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). In time domain IP surveys, we generally compute the partial chargeability (expressed in milliseconds) and defined as   \begin{equation}{M_{{t_i},{t_{i + 1}}}} = \frac{1}{{{\psi _0}}}\int_{{{t_i}}}^{{{t_{i + 1}}}}{{{\psi _{MN}}\left( t \right){\rm{d}}t}},\end{equation} (8)where $${M_{{t_i},{t_{i + 1}}}}$$is the partial chargeability measured during the time window [ti, ti + 1], ψMN(t) describes the decaying voltage measured just after the current is shut off (reference time) and ψ0 denotes the primary voltage between the potential electrodes M and N at the end of the current injection (Fig. 1a). Likewise, partial chargeability (unitless) can be expressed in mV V− 1 as:   \begin{equation}{M_{{t_i},{t_{i + 1}}}} = \frac{1}{{{\psi _0}\left( {{t_{i + 1}} - {t_i}} \right)}}\int_{{{t_i}}}^{{{t_{i + 1}}}}{{{\psi _{MN}}\left( t \right){\rm{d}}t}},\end{equation} (9)where the primary voltage is given in V and the secondary voltage in mV. Partial chargeability can be related to the apparent chargeability Ma through the approximation:   \begin{equation}{M_{{t_i},{t_{i + 1}}}} \approx {M_a}\left( {{t_{i + 1}} - {t_i}} \right).\end{equation} (10) 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 We assume now that the primary current Jp = σ0E (at low frequencies ∇ × E = 0 and therefore E = −∇ψ) has been applied from −T to time 0 (Fig. 1a). During the injection of the primary current, each cell of the discretized subsurface will see a secondary current building up. If each cell is characterized by four Cole–Cole parameters (σ0, M, τ, c), the secondary source current density is determined by (Seigel 1959),   \begin{equation}{{\bf J}_S}(t) = - M(t){{\bf J}_p}.\end{equation} (11) 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 quasi-static limit of the electromagnetic equations (taking all time derivatives to zero in the continuity equations), for each cell the constitutive equation and conservation equations are simply,   \begin{equation}{\bf J} = {{\bf J}_p} + {{\bf J}_S},\end{equation} (12)  \begin{equation}\nabla \cdot {\bf J} = I{\delta _3}({\boldsymbol{r}} - {{\boldsymbol{r}}_0}).\end{equation} (13) The fact that the total current density is solenoidal (i.e. divergence free) outside the source of primary current sources or sinks (i.e. at electrodes A and B) was extensively discussed by Seigel (1959). For each cell, the total current density can be written as,   \begin{equation}{\bf J} = \left[ {1 - M(t)} \right]{{\bf J}_p},\end{equation} (14)  \begin{equation}{\bf J} = - \left[ {1 - M(t)} \right]{\sigma _\infty }\nabla \psi .\end{equation} (15) When the primary current has been applied for long time so that the material is entirely polarized, we have M(t) = M, and we have.   \begin{equation}{\bf J} = - {\sigma _0}\nabla \psi ,\end{equation} (16)that is, the conductivity has reached its steady sate value σ0 = σ∞(1 − M). In this case, resistivity tomography provides the DC conductivity distribution of the subsurface. Our goal is to build the function M(t) during the charging phase obeying to the following properties,   \begin{equation}M(0) = 0,\end{equation} (17)  \begin{equation}M( + \infty ) = M.\end{equation} (18) For a Debye model, we have   \begin{equation}M(t) = M\left[ {1 - \exp \left( { - \frac{t}{\tau }} \right)} \right],\end{equation} (19)which is the classical behaviour for an RC circuit. For a Cole–Cole model, we have   \begin{equation}M(t) = MG\left[ {\frac{t}{\tau };c} \right],\end{equation} (20)  \begin{equation}G\left[ {\frac{t}{{{\tau _0}}};c} \right] = 1 - \sum\limits_{n = 0}^\infty {\frac{{{{( - 1)}^n}{{\left( {\frac{t}{{{\tau _0}}}} \right)}^{nc}}}}{{\Gamma (1 + nc)}}} ,\end{equation} (21)  \begin{equation}G\left[ {\frac{t}{{{\tau _0}}};c} \right] = \sum\limits_{n = 1}^\infty {\frac{{{{( - 1)}^{n + 1}}{{\left( {\frac{t}{{{\tau _0}}}} \right)}^{nc}}}}{{\Gamma (1 + nc)}}} ,\end{equation} (22)where Γ( . ) denotes the Euler gamma function defined by,   \begin{equation}\Gamma (x) = \int_{0}^{\infty }{{{u^{x - 1}}}}{e^{ - u}}{\rm{d}}u,\end{equation} (23)where x > 0. Note that the series development in eq. (22) converges very slowly for t > 10 τ and c < 1. It is easy to show that for c = 1, we recover the Debye model,   \begin{equation}G\left[ {\frac{t}{\tau };c = 1} \right] = 1 - \sum\limits_{n = 0}^\infty {\frac{{{{( - 1)}^n}{{\left( {\frac{t}{\tau }} \right)}^n}}}{{\Gamma (1 + n)}}} ,\end{equation} (24)  \begin{equation}G\left[ {\frac{t}{\tau };c = 1} \right] = 1 - \sum\limits_{n = 0}^\infty {\frac{{{{\left( { - \frac{t}{\tau }} \right)}^n}}}{{n!}}} ,\end{equation} (25)  \begin{equation}G\left[ {\frac{t}{\tau };c = 1} \right] = 1 - \exp \left( { - \frac{t}{\tau }} \right).\end{equation} (26) In the next section, we study the material response in the discharging phase. 2.3 Forward modelling in the discharging phase We first assume that the primary current has been applied from −∞ to time 0 so that the material has been entirely polarized. In the discharging phase, the primary current is zero and there is only the secondary (source) current density distribution to generate the observed electrical field, which is decaying over time. In this phase and as discussed in Mao & Revil (2016), the problem is like a transient self-potential problem given by,   \begin{equation}\nabla \cdot {\bf J} = 0,\end{equation} (27)  \begin{equation}{\bf J} = - {\sigma _\infty }\nabla \psi + {{\bf J}_S},\end{equation} (28)where for each element of the discretization characterized by the Cole–Cole parameters (σ∞, M, τ, c). We first discus the solution for a Debye model for which the source current density is simply given by   \begin{equation}{{\bf J}_S}(t) = {{\bf J}_S}(0)\exp\left( { - \frac{t}{\tau }} \right),\end{equation} (29)where the source current density at time zero is given from the end of the charging phase as,   \begin{equation}{{\bf J}_S}(0) = - M\,{{\bf J}_p},\end{equation} (30)where Jp was the primary current applied in the charging phase. Therefore, the normalized primary current in the discretized element can be written as,   \begin{equation}{{\bf \tilde{J}}_S}(t) = \frac{{{{\bf J}_S}(t)}}{{{{\bf J}_p}}} = M\exp\left( { - \frac{t}{\tau }} \right).\end{equation} (31) For the Cole–Cole model, the generalization is straightforward,   \begin{equation}{{\bf \tilde{J}}_S}(t) = M\sum\limits_{n = 0}^\infty {\frac{{{{( - 1)}^n}{{\left( {\frac{t}{\tau }} \right)}^{nc}}}}{{\Gamma (1 + nc)}}} .\end{equation} (32) If the current has been applied from time −T to zero, we can use the superposition principle by superposing a negative current from −∞ to T and a positive current of same amplitude from −∞ to zero, therefore, the solution for the discretized element   \begin{equation}{{\bf \tilde{J}}_S}(t) = M\left( {\sum\limits_{n = 0}^\infty {\frac{{{{( - 1)}^n}{{\left( {\frac{t}{\tau }} \right)}^{nc}}}}{{\Gamma (1 + nc)}}} - \sum\limits_{n = 0}^\infty {\frac{{{{( - 1)}^n}{{\left( {\frac{{t - T}}{\tau }} \right)}^{nc}}}}{{\Gamma (1 + nc)}}} } \right).\end{equation} (33) From this equation, we can compute for any element and any time the source current density distribution and then compute the resulting electrical field using,   \begin{equation}\nabla \cdot \left[ {{\sigma _\infty }\nabla \psi (t)} \right] = \nabla \cdot {{\bf J}_S}(t).\end{equation} (34) 3 COLE–COLE PARAMETERS ESTIMATION In this section, we present step-by-step the methodology that we are using for recovering the distributions of the Cole–Cole parameters by inverting the secondary voltage distribution for a set of current injection. Pelton et al. (1978) gave the Cole–Cole model formulation in time domain as   \begin{equation}M(t) = {M_0}\sum\limits_{n = 0}^\infty {\frac{{{{\left( { - 1} \right)}^n}{{\left( {\frac{t}{\tau }} \right)}^{nc}}}}{{\Gamma (1 + nc)}}} ,\end{equation} (35)where M(t) (unitless) is the intrinsic chargeability, M0(unitless) is the intrinsic chargeability at t = 0, t(s) is the time, τ(s) is the relaxation time, c(unitless) is the so-called frequency exponent and Γis the Euler Gamma function. With the current IP equipment, it is only possible to measure the apparent chargeability or secondary voltage decay curves and usually people use these curves to obtain apparent pseudo-sections of τ and c. Such apparent parameters represent averages of the intrinsic ones and may indeed be very different from it, especially if we are dealing with highly heterogeneous media. In the present work, we are looking to map in 3-D the spatial distributions of the intrinsic Cole–Cole parameters for each cell. Note that eq. (35) can have too slow convergence for large values of t/τ. Although we opted for eq. (35), alternatively, other formulations which have a faster convergence rate when t > τ can be used (e.g. Lee 1981; Hilfer 2002). The first Cole–Cole parameter that we estimate is the DC or low frequency electrical conductivity σ0. This is a classical non-linear inverse problem for which we minimize the following objective function:   \begin{eqnarray}{L_\sigma } \!=\! \left( {{\bf d}_\sigma ^{{\rm{obs}}} \!-\! {F_\sigma }\left( {\bf s} \right)} \right){\bf R}_\sigma ^{ - 1}{\left( {{\bf d}_\sigma ^{{\rm{obs}}} \!-\! {F_\sigma }\left( {\bf s} \right)} \right)^T} \!+\! \lambda \left( {{\bf s} \!-\! {{\bf s}_0}} \right){\bf C}{\left( {{\bf s} \!-\! {{\bf s}_0}} \right)^T},\end{eqnarray} (36)with $${\bf d}_\sigma ^{{\rm{obs}}}$$ denotes the (nσ × 1)observed data vector, where nσ is the number of measurements, in this case it represents the measured resistances or the apparent conductivities, Fσ(.) is the forward problem operator given by the Poisson equation, s denotes the (m × 1) model vector (unknown DC conductivities s = log 10(σ0)), and m denotes the number of unknown cells, in our case, $${\bf R}_\sigma ^{}$$is the (nσ × nσ) data covariance matrix, s0 is an a priori conductivity model, λ is the regularization parameter and can be chosen using a trial-and-error process, or some approaches such as, the L-curve approach (e.g. Hansen 1998), the Generalized Cross-validation (GCV) approach (e.g. Arlot & Celisse 2010). In our case, we will not use any prior information regarding the conductivity structure, so s0 is the null vector. The matrix C is a smoothing matrix. In our case, we chose this matrix to be a combination of first derivative operators in each space direction, that is,   \begin{equation}{\bf C} = {{\bf C}_x}{\bf C}_x^T + {{\bf C}_y}{\bf C}_y^T + {{\bf C}_z}{\bf C}_z^T,\end{equation} (37)where Cx, Cy and Cz denote the x-, y- and z-first-order derivative matrices. This could allow in future works to use an image-guided inversion to the problem (for instance for ore bodies localized along faults of know direction, see Zhou et al.2014). The next step is to use the Gauss–Newton method to minimize the cost function Lσ. With the Gauss-Newton approach, the best parameter set at a given iteration i, is updated as follows:   \begin{equation}{\bf s}_{i} = {\bf s}_{{i} - 1} + {\boldsymbol \delta}{\bf s},\end{equation} (38)with the perturbation of the model vector $${\boldsymbol \delta}{\bf s}$$ given by   \begin{equation}{\boldsymbol \delta}{\bf s} = {\left( {{{\bf J}^{\rm{T}}}{\bf R}_\sigma ^{ - 1}{\bf J}{\rm{\, +\, \lambda }}{\bf C}} \right)^{ - 1}}{{\bf J}^{\rm{T}}}{\bf R}_\sigma ^{ - 1}\left[ {{\rm{F(}}{{\boldsymbol{s}}_{{{i - 1}}}}{\rm{) - }}{\bf d}_\sigma ^{{\rm{obs}}}} \right],\end{equation} (39)and where J denotes the Jacobian matrix defined as   \begin{equation}{{\bf J}_{ij}} = \frac{{\partial {{\bf d}_i}}}{{\partial {{\bf s}_j}}}.\end{equation} (40) 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. Once a conductivity model has been obtained, we move on to the estimation of the time-dependent chargeability M(t) for each cell of the discretized subsurface, which is defined, as explained above, as,   \begin{equation}M(t) = - \frac{{{{\bf J}_S}(t)}}{{{{\bf J}_p}}},\end{equation} (41)where the primary current Jp is known and the current density distribution JS must be computed for each time in order to recover M(t) using eq. (41). We propose a novel approach for retrieving this function based on eq. (41). This inversion of the secondary current density JS can be formulated as a linear inverse problem and does not require the use of any iterative process, contrarily to the electrical conductivity problem. In fact, we seek to retrieve at each time, the amplitude and the direction of the current density distribution JS that generates the (secondary) electrical potential after the primary current has been shut-off. In order to achieve this, we discretize the simulation domain into m cells (we use the same cells as for the electrical conductivity problem) and to each cell we assign a source current density. The secondary potential decay recorded on a set of electrodes is used as data in order to invert the source current density JS. The current density at a point A of the simulation domain Ω, Js(A), can be linked to the electrical potential at a set of electrodes E, ψ(E), through:   \begin{equation}\psi \left( E \right) = \int_{\Omega }{{{\bf G}(E,A){{\bf J}_s}}}\left( A \right){\rm{d}}\Omega ,\end{equation} (42)where G(E, A) denotes the kernel matrix (a collection of Green functions), which connects the potential measured at the electrodes E to the source current density evaluated at a point A. Numerically speaking, G can be, for instance, assembled column by column, by computing the potential related to an elementary source current density at each cell (details regarding the computational aspects of the kernel can be found for instance in Trujillo-Barreto et al. 2004; Jardani et al. 2008; Soueid Ahmed et al. 2013). 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. This inverse problem is ill-posed and usually, the number of unknowns m outnumbers the number of measurements nϕ. Consequently, the inverse problem needs to be constrained to reduce the number of solutions and then to pick the optimal solution that reproduces the observed potential data and reflects the main structures of the medium as well. We point out however that the induced polarization problem we are solving is not as ill-posed as a self-potential problem. Indeed for a self-potential problem, all the source current distribution is occurring in the ground at the same time and we measure the resulting electrical potential distribution at the ground surface (plus eventually in few wells). Like all the potential field problems, this leads to a strong non-uniqueness in the solution. In the present case, the source current distribution is developed piece by piece depending on how the primary current is injected. We have therefore much more information to retrieve the distribution of the induced polarization parameters. We formulate the inverse problem as an optimization problem, for which we seek to minimize the following objective function (Tikhonov & Arsenin 1977):   \begin{equation}{L_{{J_s}}} = \left\| {{{\bf W}_d}\left( {{\bf G}{{\bf m}_{{J_S}}} - {{\bf d}^{{\rm{obs}}}}} \right)} \right\|_2^2 + \beta \big \| {{{\bf W}_m}\big( {{{\bf m}_{{J_S}}} - {{\bf m}_{J_S^0}}} \big)} \big\|_2^2,\end{equation} (43)where Wd is the diagonal nφ × nφ data covariance matrix, $${{\bf m}_{{J_S}}} = ( {{{\bf J}_{{S_x}}}\!{,}\,{{\bf J}_{{S_y}}}\!{,}\,{{\bf J}_{{S_z}}}} )$$ is the 3 m × 1 source current density model vector, dobs is the nϕ × 1vector of observed potentials, and β is the regularization parameter, Wm is a 3 m × 3 m constraint matrix and $${{\bf m}_{J_S^0}}$$ is a prior source current density model. In absence of prior information, this vector is null. However, if we have performed already some inversion with other bipoles current injection [A, B], we can formulate some prior solutions for this vector and that can further refined through the process. One should be also cautious about the choice of Wm as it plays a major role in obtaining a physically plausible solution that is clean of undesirable artefacts (see discussion in Jardani et al.2008). In our case, Wm is given by smoothness constraints defined by the x-,  y- and z-first order derivatives. The optimal solution of (43) is given by   \begin{eqnarray}{\bf m}_{{J_S}}^* &=& {\left[ {{{\bf G}^T}\left( {{\bf W}_d^T{{\bf W}_d}} \right){\bf G} + \beta \left( {{\bf W}_m^T{{\bf W}_m}} \right)} \right]^{ - 1}}\nonumber\\ &&\times \left[ {{{\bf G}^T}\left( {{\bf W}_d^T{{\bf W}_d}} \right){{\bf d}^{\rm obs}} + \beta \left( {{\bf W}_m^T{{\bf W}_m}} \right){{\bf m}_{J_S^0}}} \right].\end{eqnarray} (44) Generally, for this kind of tomography, as we get far from the sources, we lose resolution in the estimated tomograms, therefore, it is judicious to use a depth weighting matrix, which assigns weights to the cells so that, we have similar probability of having non-zeros values on the model cells, notwithstanding the fact that they have different depths (see an example of application in Haas et al. 2013). A depth weighting matrix can be computed from the kernel matrix as (e.g. Spinelli 1999; Jardani et al. 2008)   \begin{equation}{\boldsymbol{\Pi }} = {\rm{diag}}\left( {\frac{1}{{{n_\phi }}}\sqrt {\sum\limits_{j = 1}^{{n_\phi }} {{{\left( {{{\bf G}_{ij}}} \right)}^2}} } } \right)\end{equation} (45) Thus, the solution of (43) is now given by   \begin{eqnarray}{\bf m}_{J_S^\Pi }^* &=& {\left[ {{\bf G}_\Pi ^T\left( {{\bf W}_d^T{{\bf W}_d}} \right){{\bf G}_\Pi } + \beta \left( {{\bf W}_m^T{{\bf W}_m}} \right)} \right]^{ - 1}}\nonumber\\ &&\times \left[ {{\bf G}_\Pi ^T\left( {{\bf W}_d^T{{\bf W}_d}} \right){{\bf d}^{{\rm{obs}}}} + \beta \left( {{\bf W}_m^T{{\bf W}_m}} \right){{\bf m}_{J_S^0}}} \right]\end{eqnarray} (46)where   \begin{equation}{{\bf G}_\Pi } = {\bf G}{{\boldsymbol \Pi }^{ - 1}}.\end{equation} (47) We retrieve $${\bf m}_{{J_S}}^*$$, by writing   \begin{equation}{\bf m}_{{J_S}}^* = {{\bf m}_{{J_S}}}{{\boldsymbol{\Pi }}^{ - 1}}\end{equation} (48) 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. Figure 2. View largeDownload slide Step-by-step summary of the inversion methodology. ERT corresponds to electrical resistivity tomography. The first step is to determine the electrical conductivity field from which the kernel used to invert the secondary voltages is computed. For each injected current at a bipole (A, B), we compute the secondary source current density using a linear inversion procedure (second step). The chargeability is given as the ratio between the secondary and the primary current distribution. Finally the relaxation time and the Cole–Cole frequency exponent are determined at each cell by fitting the inverted data and a variety of strategies can be used for this final step (third step). Figure 2. View largeDownload slide Step-by-step summary of the inversion methodology. ERT corresponds to electrical resistivity tomography. The first step is to determine the electrical conductivity field from which the kernel used to invert the secondary voltages is computed. For each injected current at a bipole (A, B), we compute the secondary source current density using a linear inversion procedure (second step). The chargeability is given as the ratio between the secondary and the primary current distribution. Finally the relaxation time and the Cole–Cole frequency exponent are determined at each cell by fitting the inverted data and a variety of strategies can be used for this final step (third step). 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). Figure 3. View largeDownload slide Geometry of the synthetic experiment. For this case, we used a 0.5m × 0.5m × 0.5m domain with insulating boundary conditions (simulating a sandbox experiment). In this domain, we place 24 electrodes denoted by the dots in the figure. These electrodes are placed in such a way that they encircle the central arc of the tank, where the polarizable body is placed. For the purpose of this experiment, 2 current injection/retrieval pair [A, B] have been performed at the bipoles [23, 28] and [7, 10], respectively. The electrode #1 is used as reference electrode for the secondary voltages. This box is discretized with 1000 elements. Figure 3. View largeDownload slide Geometry of the synthetic experiment. For this case, we used a 0.5m × 0.5m × 0.5m domain with insulating boundary conditions (simulating a sandbox experiment). In this domain, we place 24 electrodes denoted by the dots in the figure. These electrodes are placed in such a way that they encircle the central arc of the tank, where the polarizable body is placed. For the purpose of this experiment, 2 current injection/retrieval pair [A, B] have been performed at the bipoles [23, 28] and [7, 10], respectively. The electrode #1 is used as reference electrode for the secondary voltages. This box is discretized with 1000 elements. Figure 4. View largeDownload slide True and estimated Cole–Cole parameters distributions. (a) True Cole–Cole parameters models. (b) Estimated electrical conductivity at iteration 4 (RMS = 0.8, convergence criteria: variation of the objective function < 0.001 between two successive iterations). A 5 per cent Gaussian noise has been added to the synthetic data. (c) Estimated intrinsic chargeability at time 0. (d) Estimated relaxation time. (e) Estimated frequency exponent. We can note that the anomaly is well recovered both in terms of magnitude and location. Figure 4. View largeDownload slide True and estimated Cole–Cole parameters distributions. (a) True Cole–Cole parameters models. (b) Estimated electrical conductivity at iteration 4 (RMS = 0.8, convergence criteria: variation of the objective function < 0.001 between two successive iterations). A 5 per cent Gaussian noise has been added to the synthetic data. (c) Estimated intrinsic chargeability at time 0. (d) Estimated relaxation time. (e) Estimated frequency exponent. We can note that the anomaly is well recovered both in terms of magnitude and location. Table 1. Synthetic IP experiment data configuration. In the synthetic experiment, we inject a current of 1 mA during 2 s and we record the secondary voltage decay after a dead time of 0.1 s. 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  View Large Table 2. Sandbox experiment IP data configuration. The current injected between two electrodes A and B, for 1 s and the secondary voltages are recorded on all the remaining electrodes, taking electrode #1 as a reference. Two current injections are performed during this experiment. After the current is shut down, we wait 0.03 s before starting recording the secondary voltages. This delay time is used to avoid the electromagnetic effects on IP data. The last time, that is, 5.13 s is used as a temporal reference for all the measured voltages. 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.  View Large 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. Figure 5. View largeDownload slide Primary and secondary current density distributions. (a) Primary current density when the current electrodes A and B correspond to the red electrodes shown in Fig. 3. The current came across the target and its density quickly away from A and B. (b) Secondary (source) current density field. This field reflects an area of high current density, which corresponds to the target. Figure 5. View largeDownload slide Primary and secondary current density distributions. (a) Primary current density when the current electrodes A and B correspond to the red electrodes shown in Fig. 3. The current came across the target and its density quickly away from A and B. (b) Secondary (source) current density field. This field reflects an area of high current density, which corresponds to the target. Figure 6. View largeDownload slide Intrinsic chargeability decay during each time window after the primary current was shut off. Note that the chargeability at time 0 is obtained jointly with the relaxation time and the frequency exponent using the curve fitting method. It cannot be inferred from the linear inversion because the potential is not recorded at time 0 since early times are contaminated with electromagnetic coupling effects. Figure 6. View largeDownload slide Intrinsic chargeability decay during each time window after the primary current was shut off. Note that the chargeability at time 0 is obtained jointly with the relaxation time and the frequency exponent using the curve fitting method. It cannot be inferred from the linear inversion because the potential is not recorded at time 0 since early times are contaminated with electromagnetic coupling effects. 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). Figure 7. View largeDownload slide Sketch of the sandbox. The tank is filled with water-saturated sand and contains a target composed of burning coal placed in the cylindrical area located at the central part of the sandbox. The dots denote the 52 stainless electrodes that are used for data acquisition. Two current injections are performed at the electrodes pairs [A, B] = [5, 14] (red) and [31, 40] (blue), respectively. The electrode #1 is used as the reference electrode (Ref) for the voltage electrodes recording the secondary voltage decay. The tank is discretized with 1000 elements. Figure 7. View largeDownload slide Sketch of the sandbox. The tank is filled with water-saturated sand and contains a target composed of burning coal placed in the cylindrical area located at the central part of the sandbox. The dots denote the 52 stainless electrodes that are used for data acquisition. Two current injections are performed at the electrodes pairs [A, B] = [5, 14] (red) and [31, 40] (blue), respectively. The electrode #1 is used as the reference electrode (Ref) for the voltage electrodes recording the secondary voltage decay. The tank is discretized with 1000 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). Figure 8. View largeDownload slide Relaxation of the secondary voltages following the shutdown of the primary current associated with the bipole [A, B] = [5 14]. (a) Secondary voltage recorded at electrodes #2 to #9 (see position in Fig. 7). (b) Secondary voltage recorded at electrodes #9 to #15 (see position in Fig. 7). The secondary voltages exhibit some strong IP effects, which are associated with the presence of the coal (the same experiment performed without coal does not exhibit such polarization effect). The voltage at t5 = 13 s is used as a temporal reference for the secondary voltages at the different measurement times (i.e. we consider that the secondary voltages has fully relaxed at this time). Figure 8. View largeDownload slide Relaxation of the secondary voltages following the shutdown of the primary current associated with the bipole [A, B] = [5 14]. (a) Secondary voltage recorded at electrodes #2 to #9 (see position in Fig. 7). (b) Secondary voltage recorded at electrodes #9 to #15 (see position in Fig. 7). The secondary voltages exhibit some strong IP effects, which are associated with the presence of the coal (the same experiment performed without coal does not exhibit such polarization effect). The voltage at t5 = 13 s is used as a temporal reference for the secondary voltages at the different measurement times (i.e. we consider that the secondary voltages has fully relaxed at this time). 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. Figure 9. View largeDownload slide Tomogram of the electrical conductivity field at the fourth iteration (RMS = 5.7, convergence criteria: variation of the objective function <0.001 between two successive iterations). The tomogram shows an area of high conductivity which is located exactly where the coal is burning. This is consistent with the fact that burning coal exhibits a strong conductivity. Figure 9. View largeDownload slide Tomogram of the electrical conductivity field at the fourth iteration (RMS = 5.7, convergence criteria: variation of the objective function <0.001 between two successive iterations). The tomogram shows an area of high conductivity which is located exactly where the coal is burning. This is consistent with the fact that burning coal exhibits a strong conductivity. Figure 10. View largeDownload slide Observed versus computed resistances at the position of the 50 electrodes. (a) Observed versus computed resistances for bipole [A, B] = [31, 40]. (b) Observed versus computed resistances for bipole [A, B] = [5, 14]. The computed resistances reproduce the measured resistances. The latter were used as data for the electrical resistivity tomography. Figure 10. View largeDownload slide Observed versus computed resistances at the position of the 50 electrodes. (a) Observed versus computed resistances for bipole [A, B] = [31, 40]. (b) Observed versus computed resistances for bipole [A, B] = [5, 14]. The computed resistances reproduce the measured resistances. The latter were used as data for the electrical resistivity tomography. Figure 11. View largeDownload slide Primary current density distributions. (a) Primary current distribution for the electrode bipole [A, B] = [31; 40]. (b) Primary current distribution for the electrode bipole [A, B] = [5 14]. The primary current distributions are very high in the vicinity of the sources, that is, around the current injection electrodes. We observe the presence pattern of sensitivity of the primary current density, and as we get far from this pattern, the sensitivity drops and the current density becomes very small. Figure 11. View largeDownload slide Primary current density distributions. (a) Primary current distribution for the electrode bipole [A, B] = [31; 40]. (b) Primary current distribution for the electrode bipole [A, B] = [5 14]. The primary current distributions are very high in the vicinity of the sources, that is, around the current injection electrodes. We observe the presence pattern of sensitivity of the primary current density, and as we get far from this pattern, the sensitivity drops and the current density becomes very small. We recover now the source current density field on each time window, from the secondary voltage potentials. As stated before in Section 2, this is a linear inverse problem which can be solved without any iterative process. The kernel matrix which relates the source current density field to the voltage data is computed once and then used to retrieve the optimal source current density at each measurement time using eq. (46). As an example, we represent on Fig. 12 the source current density field obtained at time 0.16 s. It reveals that we have a high current density anomaly located at the central region of the tank. This indicates the presence of a highly polarizable body, which indeed coincides with the area where we ignited the coal. Fig. 13 shows that the estimated source current density reproduces with high fidelity the observed voltages. The next step is to recover the intrinsic chargeability M(t) distribution. We estimate M(t), where ‖ . ‖denotes the norm operator. In order to achieve this, we write   \begin{equation}\left\| {M\left( t \right)} \right\| = \frac{{\left\| {{{\bf J}_S}\left( t \right)} \right\|}}{{\Vert {{{\bf J}_p}}\Vert}}.\end{equation} (49) Figure 12. View largeDownload slide Secondary current density distribution. This distribution is here estimated at t = 0.16 s. It shows an anomaly of high current density which indeed is located in the central area where the burning coal is located. Figure 12. View largeDownload slide Secondary current density distribution. This distribution is here estimated at t = 0.16 s. It shows an anomaly of high current density which indeed is located in the central area where the burning coal is located. Figure 13. View largeDownload slide Observed versus computed secondary voltages. (a) Comparison for bipole [A, B] = [31, 40]. (b) Comparison for bipole [A, B] = [5, 14]. We plot the true and estimated secondary voltages on all the electrodes at time t = 0.16 s. The computed voltages are obtained using the estimated source current density corresponding to the optimal solution of linear inverse problem. Figure 13. View largeDownload slide Observed versus computed secondary voltages. (a) Comparison for bipole [A, B] = [31, 40]. (b) Comparison for bipole [A, B] = [5, 14]. We plot the true and estimated secondary voltages on all the electrodes at time t = 0.16 s. The computed voltages are obtained using the estimated source current density corresponding to the optimal solution of linear inverse problem. 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. Figure 14. View largeDownload slide Intrinsic chargeability as a function of the elapsed time after the shutdown of the primary current. (a–h) We observe the intrinsic chargeability temporal evolution just after the primary current is shut down. The intrinsic chargeability exhibits a high anomaly in the central area of the tank, which is associated with the burning coal. The amplitude of the anomaly keeps decreasing over time, until it completely vanishes. Figure 14. View largeDownload slide Intrinsic chargeability as a function of the elapsed time after the shutdown of the primary current. (a–h) We observe the intrinsic chargeability temporal evolution just after the primary current is shut down. The intrinsic chargeability exhibits a high anomaly in the central area of the tank, which is associated with the burning coal. The amplitude of the anomaly keeps decreasing over time, until it completely vanishes. 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). Figure 15. View largeDownload slide Initial chargeability, relaxation time and Cole–Cole exponent fields. (a) Initial chargeability. (b) Relaxation time (c) Frequency exponent. We represent τ and c only in the region characterized by high values of M0 (>0.15). Indeed, the distribution of M0 shows high polarization effects where the coal is burning. The τ values are the highest (around 1.5 s) in the region where the coal is burning, while the values of c remain relatively high (greater than 0.9), suggesting that the polarization of the burning coal can be approximated with a Debye model. Figure 15. View largeDownload slide Initial chargeability, relaxation time and Cole–Cole exponent fields. (a) Initial chargeability. (b) Relaxation time (c) Frequency exponent. We represent τ and c only in the region characterized by high values of M0 (>0.15). Indeed, the distribution of M0 shows high polarization effects where the coal is burning. The τ values are the highest (around 1.5 s) in the region where the coal is burning, while the values of c remain relatively high (greater than 0.9), suggesting that the polarization of the burning coal can be approximated with a Debye model. 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. Figure 16. View largeDownload slide Inverted chargeability versus time in two arbitrary cells belonging to the computation domain. (a) Cell 1. (b) Cell 2. We chose these two arbitrary cells of the simulation domain in which we plot the intrinsic chargeability decay (blue dots) and its fitting with a time domain Cole–Cole relaxation model. By fitting the chargeability decay in each cell, we recover a value of M0, τand con the corresponding cell. Figure 16. View largeDownload slide Inverted chargeability versus time in two arbitrary cells belonging to the computation domain. (a) Cell 1. (b) Cell 2. We chose these two arbitrary cells of the simulation domain in which we plot the intrinsic chargeability decay (blue dots) and its fitting with a time domain Cole–Cole relaxation model. By fitting the chargeability decay in each cell, we recover a value of M0, τand con the corresponding cell. 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 Arlot A., Celisse A., 2010. A survey of cross-validation procedures for model selection, Stat. Surv.  4( 0), 40– 79. https://doi.org/10.1214/09-SS054 Google Scholar CrossRef Search ADS   Bérubé C.L., Chouteau M., Shamsipour P., Enkin R.J., Olivo G.R., 2017. Bayesian inference of spectral induced polarization parameters for laboratory complex resistivity measurements of rocks and soils, Comput. Geosci.  105 51– 64. https://doi.org/10.1016/j.cageo.2017.05.001 Google Scholar CrossRef Search ADS   Binley A., Hubbard S.S., Huisman J.A., Revil A., Robinson D.A., Singha K., Slater L.D., 2015. The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales, Water Resour. Res.  51( 6), 3837– 3866. https://doi.org/10.1002/2015WR017016 Google Scholar CrossRef Search ADS PubMed  Bleil D.F., 1953, Induced polarization: A method of geophysical prospecting, Geophysics  18( 3), 636– 661. https://doi.org/10.1190/1.1437917 Google Scholar CrossRef Search ADS   Chen J., Kemna A., Hubbard S.S., 2008. A comparison between Gauss–Newton and Markov-chain Monte Carlo–based methods for inverting spectral induced-polarization data for Cole–Cole parameters, Geophysics  73( 6), F247– F259. https://doi.org/10.1190/1.2976115 Google Scholar CrossRef Search ADS   Cole K.S., Cole R.H., 1941. Dispersion and absorption in dielectrics I. Alternating current characteristics, J. Chem. Phys.  9( 4), 341– 351. https://doi.org/10.1063/1.1750906 Google Scholar CrossRef Search ADS   Comte J.C., Banton O., Join J.-L., Cabioch G., 2010. Evaluation of effective groundwater recharge of freshwater lens in small islands by the combined modeling of geoelectrical data and water heads, Water Resour. Res.  46( 6), W06601, doi:10.1029/2009WR008058. https://doi.org/10.1029/2009WR008058 Google Scholar CrossRef Search ADS   Dakhnov V.N., 1962, Geophysical well logging, Colorado School Mines Quart.  57( 2), 445. Fiandaca G., Auken E., Christiansen A.V., Gazoty A., 2012. Time-domain-induced polarization: full-decay forward modeling and 1D laterally constrained inversion of Cole-Cole parameters, Geophysics  77( 3), E213– E225. https://doi.org/10.1190/geo2011-0217.1 Google Scholar CrossRef Search ADS   Fiandaca G., Ramm J., Binley A., Gazoty A., Christiansen A.V., Auken E., 2013. Resolving spectral information from time domain induced polarization data through 2-D inversion, Geophys. J. Int.  192( 2), 631– 646. https://doi.org/10.1093/gji/ggs060 Google Scholar CrossRef Search ADS   Florsch N., Camerlynck C., Revil A., 2012. Direct estimation of the distribution of relaxation times from induced-polarization spectra using a Fourier transform analysis. Near Surface Geophysics  10( 6), 517– 531. https://doi.org/10.1093/gji/ggs060 Florsch N., Llubes M., Téreygeol F., Ghorbani A., Roblet P., 2011. Quantification of slag heap volumes and masses through the use of induced polarization: application to the Castel-Minier site, J. Archaeological Sci.  38( 2), 438– 451. https://doi.org/10.1016/j.jas.2010.09.027 Google Scholar CrossRef Search ADS   Ghorbani A., Camerlynck C., Florsch N., Cosenza P., Revil A. 2007. Bayesian inference of the Cole–Cole parameters from time- and frequency-domain induced polarization, Geophys. Prospect.  55( 4), 589– 605. https://doi.org/10.1111/j.1365-2478.2007.00627.x Google Scholar CrossRef Search ADS   Haas A., Revil A., Karaoulis M., Frash L., Hampton J., Gutierrez M., Mooney M., 2013. Electric potential source localization reveals a borehole leak during hydraulic fracturing, Geophysics  78( 2), D93– D113. https://doi.org/10.1190/geo2012-0388.1 Google Scholar CrossRef Search ADS   Hansen P.C., 1998. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion  SIAM, p. 247. Google Scholar CrossRef Search ADS   Hilfer R., 2002. Analytical representations for relaxation functions of glasses, J. Non-Cryst. Solids  305( 1-3), 122– 126. https://doi.org/10.1016/S0022-3093(02)01088-8 Google Scholar CrossRef Search ADS   Ikard S.J., Revil A., Jardani A., Woodruff W.F., Parekh M., Mooney M., 2012. Saline pulse test monitoring with the self-potential method to nonintrusively determine the velocity of the pore water in leaking areas of earth dams and embankments, Water Resour. Res.  48( 4), W04201, doi:10.1029/2010WR010247. https://doi.org/10.1029/2010WR010247 Google Scholar CrossRef Search ADS   Jardani A., Revil A., Bolève A., Dupont J.P., 2008. Three-dimensional inversion of self-potential data used to constrain the pattern of groundwater flow in geothermal fields, J. geophys. Res.  113( B9), B09204, doi:10.1029/2007JB005302. https://doi.org/10.1029/2007JB005302 Google Scholar CrossRef Search ADS   Jardani A., Revil A., Dupont J.P., 2013. Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging, Adv. Water Resour.  52 62– 77. https://doi.org/10.1016/j.advwatres.2012.08.005 Google Scholar CrossRef Search ADS   Joseph S., Ingham M., Gouws G., 2016. Spectral-induced polarization measurements on sieved sands and the relationship to permeability, Water Resour. Res.  52( 6), 4226– 4246. https://doi.org/10.1002/2015WR018087 Google Scholar CrossRef Search ADS   Keery J., Binley A., Elshenawy A., Clifford J., 2012. Markov-chain Monte Carlo estimation of distributed Debye relaxations in spectral induced polarization, Geophysics  77( 2), E159– E170. https://doi.org/10.1190/geo2011-0244.1 Google Scholar CrossRef Search ADS   Kemna A., Binley A., Slater L., 2004. Crosshole IP imaging for engineering and environmental applications, Geophysics  69( 1), 97– 107. https://doi.org/10.1190/1.1649379 Google Scholar CrossRef Search ADS   Kemna A. et al.  , 2012. An overview of the spectral induced polarization method for near-surface applications, Near Surf. Geophys.  10( 1957), 453– 468. https://doi.org/10.3997/1873-0604.2012027 Lee T., 1981. The Cole–Cole model in time domain induced polarization, Geophysics  46( 6), 932– 933. https://doi.org/10.1190/1.1441231 Google Scholar CrossRef Search ADS   Loke M.H., Chambers J.E., Ogilvy R.D., 2006. Inversion of 2D spectral induced polarization imaging data, Geophys. Prospect.  54( 3), 287– 301. https://doi.org/10.1111/j.1365-2478.2006.00537.x Google Scholar CrossRef Search ADS   Mahardika H., Revil A., Jardani A., 2012. Waveform joint inversion of seismograms and electrograms for moment tensor characterization of fracking events, Geophysics  77( 5), ID23– ID39. https://doi.org/10.1190/geo2012-0019.1 Google Scholar CrossRef Search ADS   Mao D., Revil A., 2016. Induced polarization response of porous media with metallic particles—Part 3: A new approach to time-domain induced polarization tomography, Geophysics  81( 4), D345– D357. https://doi.org/10.1190/geo2015-0283.1 Google Scholar CrossRef Search ADS   Mao D., Revil A., Hinton J., 2016. Induced polarization response of porous media with metallic particles. – Part 4: Detection of metallic and nonmetallic targets in time-domain induced polarization tomography, Geophysics  81( 4), D359– D375. https://doi.org/10.1190/geo2015-0480.1 Google Scholar CrossRef Search ADS   Michot D., Benderitter Y., Dorigny A., Nicoullaud B., King D., Tabbagh A., 2003. Spatial and temporal monitoring of soil water content with an irrigated corn crop cover using surface electrical resistivity tomography, Water Resour. Res.  39( 5), 1138, doi:10.1029/2002WR001518. https://doi.org/10.1029/2002WR001581 Google Scholar CrossRef Search ADS   Misra S., Torres-Verdín C., Revil A., Rasmus J., Homan D., 2016a. Interfacial polarization of disseminated conductive minerals in absence of redox-active species—Part 1: Mechanistic model and validation, Geophysics  81( 2), E139– E157. https://doi.org/10.1190/geo2015-0346.1 Google Scholar CrossRef Search ADS   Misra S., Torres-Verdín C., Revil A., Rasmus J., Homan D., 2016b. Interfacial polarization of disseminated conductive minerals in absence of redox-active species—Part 2: Effective electrical conductivity and dielectric permittivity, Geophysics  81( 2), E159– E176. https://doi.org/10.1190/geo2015-0400.1 Google Scholar CrossRef Search ADS   Niu Q., Revil A., 2016. Connecting complex conductivity spectra to mercury porosimetry of sedimentary rocks, Geophysics  81( 1), E17– E32. https://doi.org/10.1190/geo2015-0072.1 Google Scholar CrossRef Search ADS   Nivorlis A., Tsourlos P., Vargemezis G., Tsokas G., Kim J.H., Yi M.J., 2017. Processing and Modeling of Time Domain Induced Polarization Data, 23rd European Meeting of Environmental and Engineering Geophysics, extended abstract, doi: 10.3997/2214-4609.201702085. Olhoeft G.R., 1984, Clay-organic interactions measured with complex resistivity, in 54th Annual International Meeting and Exposition of the Society of Exploration Geophysicists  Expanded Abstracts, 356– 358. Olhoeft G.R., 1985. Low?frequency electrical properties, Geophysics  50( 12), 2492– 2503. https://doi.org/10.1190/1.1441880 Google Scholar CrossRef Search ADS   Olhoeft G.R., 1986. Direct detection of hydrocarbon and organic chemicals with ground-penetrating radar and complex resistivity, in Proceedings of the NWWA/API Conference on Petroleum Hydrocarbons and Organic Chemicals in GroundWater-Prevention, Detection, and Restoration , pp. 284– 305, 1985 November 13–15, Houston, TX. Osterman G., Keating K., Binley A., Slater L., 2016. A laboratory study to estimate pore geometric parameters of sandstones using complex conductivity and nuclear magnetic resonance for permeability prediction, Water Resour. Res.  52( 6), 4321– 4337. https://doi.org/10.1002/2015WR018472 Google Scholar CrossRef Search ADS   Pelton W.H., Ward S.H., Hallof P.G., Sill W.R., Nelson P.H., 1978, Mineral discrimination and removal of inductive coupling with multifrequency IP, Geophysics  43( 3), 588– 609. https://doi.org/10.1190/1.1440839 Google Scholar CrossRef Search ADS   Revil A., 2012. Spectral induced polarization of shaly sands: Influence of the electrical double layer, Water Resour. Res.  48( 2), W02517, doi:10.1029/2011WR011260. https://doi.org/10.1029/2011WR011260 Google Scholar CrossRef Search ADS   Revil A., 2013. Effective conductivity and permittivity of unsaturated porous materials in the frequency range 1 mHz-1 GHz, Water Resour. Res.  49( 1), 306– 327. https://doi.org/10.1029/2012WR012700 Google Scholar CrossRef Search ADS PubMed  Revil A., 2017a. Transport of water and ions in partially water-saturated porous media—Part 1: Constitutive equations, Advances in Water Resources  103 119– 138. https://doi.org/10.1016/j.advwatres.2016.02.006 Google Scholar CrossRef Search ADS   Revil A., 2017b. Transport of water and ions in partially water-saturated porous media—Part 2: Filtration effects, Adv. Water Resour.  103 139– 152. https://doi.org/10.1016/j.advwatres.2016.07.016 Google Scholar CrossRef Search ADS   Revil A., Cathles L.M., Losh S., Nunn J.A., 1998. Electrical conductivity in shaly sands with geophysical applications, J. geophys. Res.  103( B10), 23 925–23 936. https://doi.org/10.1029/98JB02125 Revil A., Florsch N., Camerlynck C., 2014. Spectral induced polarization porosimetry, Geophys. J. Int.  198( 2), 1016– 1033. https://doi.org/10.1093/gji/ggu180 Google Scholar CrossRef Search ADS   Revil A., Florsch N., Mao D., 2015a. Induced polarization response of porous media with metallic particles—Part 1: A theory for disseminated semiconductors, Geophysics  80( 5), D525– D538. https://doi.org/10.1190/geo2014-0577.1 Google Scholar CrossRef Search ADS   Revil A., Abdel Aal G.Z., Atekwana E.A., Mao D., Florsch N., 2015b. Induced polarization response of porous media with metallic particles—Part 2: Comparison with a broad database of experimental data, Geophysics  80( 5), D539– D552. https://doi.org/10.1190/geo2014-0578.1 Google Scholar CrossRef Search ADS   Revil A., Binley A., Mejus L., Kessouri P., 2015c. Predicting permeability from the characteristic relaxation time and intrinsic formation factor of complex conductivity spectra, Water Resour. Res.  51( 8), 6672– 6700. https://doi.org/10.1002/2015WR017074 Google Scholar CrossRef Search ADS   Revil A., Le Breton M., Niu Q., Wallin E., Haskins E., Thomas D.M., 2017a. Induced polarization of volcanic rocks—1. Surface versus quadrature conductivity, Geophys. J. Int.  208( 2), 826– 844. https://doi.org/10.1093/gji/ggw444 Google Scholar CrossRef Search ADS   Revil A., Le Breton M., Niu Q., Wallin E., Haskins E., Thomas D.M., 2017b. Induced polarization of volcanic rocks—2. Influence of pore size and permeability, Geophys. J. Int.  208( 2), 814– 825. https://doi.org/10.1093/gji/ggw382 Google Scholar CrossRef Search ADS   Rosen L.A., Saville D.A., 1991. Dielectric spectroscopy of colloidal dispersions: comparisons between experiment and theory, Langmuir  7( 1), 36– 42. https://doi.org/10.1021/la00049a009 Google Scholar CrossRef Search ADS   Rosen L.A., Baygents J.C., Saville D.A., 1993. The interpretation of dielectric response measurements on colloidal dispersions using the dynamic Stern layer model, J. Chem. Phys.  98( 5), 4183– 4194. https://doi.org/10.1063/1.465108 Google Scholar CrossRef Search ADS   Schlumberger C., 1920, Etude Sur La Prospection Électrique Du Sous-Sol  Rev. edn, Gauthier Villars. Seigel H.O., 1949. Theoretical and experimental investigations into the applications of the phenomenon of overvoltage to geophysical prospecting, Toronto, Doctoral dissertation  University of Toronto. Seigel H.O., 1959, Mathematical formulation and type curves for induced polarization, Geophysics  24( 3), 547– 565. https://doi.org/10.1190/1.1438625 Google Scholar CrossRef Search ADS   Shainberg I., Rhoades J.D., Prather R.J., 1980, Effect of exchangeable sodium percentage, cation exchange capacity, and soil solution concentration on soil electrical conductivity1, Soil Sci. Soc. Am. J.  44( 3), 469– 473. https://doi.org/10.2136/sssaj1980.03615995004400030006x Google Scholar CrossRef Search ADS   Shao Z., Revil A., Mao D., Wang D. 2017. Induced polarization signature of coal seam fires, Geophys. J. Int.  208( 3), 1313– 1331. https://doi.org/10.1093/gji/ggw452 Google Scholar CrossRef Search ADS   Slater L.D., Lesmes D. 2002. IP interpretation in environmental investigations, Geophysics  67( 1), 77– 88. https://doi.org/10.1190/1.1451353 Google Scholar CrossRef Search ADS   Soueid Ahmed A., Jardani A., Revil A., Dupont J.P. 2013. SP2DINV: A 2D forward and inverse code for streaming potential problems, Comput. Geosci.  59 9– 16. https://doi.org/10.1016/j.cageo.2013.05.008 Google Scholar CrossRef Search ADS   Spinelli L. 1999. Analyse spatiale de l'activité électrique cérébrale: nouveaux développements, Doctoral dissertation  Université Joseph-Fourier-Grenoble I, France. Sturrock J.T., Lesmes D., Morgan F.D., 1999. Permeability estimation using spectral induced polarization measurements, in Proc. Symp. on the Application of Geophysics to Engineering and Environmental Problems , pp. 409– 415. Tarasov K., Titov K., 2013, On the use of the Cole–Cole equations in spectral induced polarization, Geophys. J. Int.  195( 1), 352– 356. https://doi.org/10.1093/gji/ggt251 Google Scholar CrossRef Search ADS   Telford W.M., Geldart L.P., Sheriff R.E., 1990. Applied Geophysics  second revised edn, p 792, Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Tikhonov A.N., Arsenin V.Y., 1977. Solutions of Ill-posed Problems  John Wiley & Sons. Trujillo-Barreto N.J., Aubert-Vásquez E., Valdès-Sosa P.A., 2004. Bayesian model averaging in EEG-MEG imaging, Neuro Image  21 1300– 1319. Google Scholar PubMed  Vanhala H., 1997, Mapping oil-contaminated sand and till with the spectral induced polarization (SIP) method, Geophys. Prospect.  45( 2), 303– 326. https://doi.org/10.1046/j.1365-2478.1997.00338.x Google Scholar CrossRef Search ADS   Vanhala H., Soininen H., Kukkonen I., 1992. Detecting organic chemical contaminants by spectral-induced polarization method in glacial till environment, Geophysics  57( 8), 1014– 1017. https://doi.org/10.1190/1.1443312 Google Scholar CrossRef Search ADS   Van Voorhis G.D., Nelson P.H., Drake T.L., 1973. Complex resistivity spectra of porphyry copper mineralization, Geophysics  38( 1), 49– 60. https://doi.org/10.1190/1.1440333 Google Scholar CrossRef Search ADS   Wainwright H.M., Chen J., Sassen D.S., Hubbard S.S., 2014, Bayesian hierarchical approach and geophysical data sets for estimation of reactive facies over plume scales, Water Resour. Res.  50( 6), 4564– 4584. https://doi.org/10.1002/2013WR013842 Google Scholar CrossRef Search ADS   Waxman M.H., Smits L.J.M., 1968. Electrical conductivities in oil-bearing shaly sands, SPE J.  8( 02), 107– 122. https://doi.org/10.2118/1863-A Google Scholar CrossRef Search ADS   Weller A., Breede K., Slater L., Nordsiek S., 2011. Effect of changing water salinity on complex conductivity spectra of sandstones, Geophysics  76( 5), F315– F327. https://doi.org/10.1190/geo2011-0072.1 Google Scholar CrossRef Search ADS   Weller A., Slater L., Nordsiek S., 2013. On the relationship between induced polarization and surface conductivity: implications for petrophysical interpretation of electrical measurements, Geophysics  78( 5), D315– D325. https://doi.org/10.1190/geo2013-0076.1 Google Scholar CrossRef Search ADS   Weller A., Slater L., Huisman J.A., Esser O., Haegel F.H., 2015. On the specific polarizability of sands and sand-clay mixtures, Geophysics  80( 3), A57– A61. https://doi.org/10.1190/geo2014-0509.1 Google Scholar CrossRef Search ADS   Weller A., Slater L., Binley A., Nordsiek S., Xu S., 2015. Permeability prediction based on induced polarization: Insights from measurements on sandstone and unconsolidated samples spanning a wide permeability range, Geophysics  80( 2), D161– D173. https://doi.org/10.1190/geo2014-0368.1 Google Scholar CrossRef Search ADS   Wong J., 1979. An electrochemical model of the induced polarization phenomenon in disseminated sulfide ores, Geophysics  44 1245– 1265. https://doi.org/10.1190/1.1441005 Google Scholar CrossRef Search ADS   Yuval, Oldenburg D.W., 1997. Computation of Cole–Cole parameters from IP data, Geophysics  62( 2), 436– 448. https://doi.org/10.1190/1.1444154 Google Scholar CrossRef Search ADS   Zimmermann E., Kemna A., Berwix J., Glaas W., Munch H., Huisman J., 2008. A high-accuracy impedance spectrometer for measuring sediments with low polarizability, Meas. Sci. Technol.  19( 10), 105603. https://doi.org/10.1088/0957-0233/19/10/105603 Google Scholar CrossRef Search ADS   Zhou J., Revil A., Karaoulis M., Hale D., Doetsch J., Cuttler S., 2014. Image-guided inversion of electrical resistivity data, Geophys. J. Int.  197( 1), 292– 309. https://doi.org/10.1093/gji/ggu001 Google Scholar CrossRef Search ADS   Zonge K.L., Wynn J.C., 1975. Recent advances and applications in complex resistivity measurements, Geophysics  40( 5), 851– 864. https://doi.org/10.1190/1.1440572 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - 3-D time-domain induced polarization tomography: a new approach based on a source current density formulation
JF - Geophysical Journal International
DO - 10.1093/gji/ggx547
DA - 2018-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/3-d-time-domain-induced-polarization-tomography-a-new-approach-based-yGIvV1IUp4
SP - 244
EP - 260
VL - 213
IS - 1
DP - DeepDyve
ER -