3-D time-domain induced polarization tomography: a new approach based on a source current density formulation

About: This article is published in Geophysical Journal International.The article was published on 2017-12-18 and is currently open access. It has received 9 citations till now. The article focuses on the topics: Induced polarization & Time domain.

Summary

Introduction

• Received 2017 December 13; in original form 2017 April 17 S U M M A R Y 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 current density.
• IP effects were first discovered during the last century (e.g. Schlumberger 1920; Dakhnov 1962).
• The advantage of formulating the IP problem as an equivalent self-potential problem is that the authors avoid solving a nonlinear inverse problem because retrieving the current density JS(t) from the recorded electrical field is a linear problem.
• To the best of their knowledge, this is the first attempt to recover the Cole–Cole parameters in 3-D using such an approach.

2.1 Electrical conductivity and chargeability

• The authors consider below time scale and length scales such as the induction effect can be neglected.
• This potential ψσ∞ is instantaneously recorded when the current injection is turned on.
• (6) Similarly, an apparent chargeability can be computed as Ma = ψσ∞ − ψσ0 ψσ∞ = F (σ∞) − F (σ∞ (1 − M)) F (σ∞) . (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).
• Note that their notations are pretty standard.

2.2 Forward modelling in the charging phase

• 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), JS(t) = −M(t)Jp. (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.
• 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, J = Jp + JS, (12) ∇ · J = I δ3(r − r0). (13).
• When the primary current has been applied for long time so that the material is entirely polarized, the authors have M(t) = M, and they have.

2.3 Forward modelling in the discharging phase

• The authors first assume that the primary current has been applied from −∞ to time 0 so that the material has been entirely polarized.
• Alternatively, other formulations which have a faster convergence rate when t > τ can be used (e.g. Lee 1981; Hilfer 2002).the authors.
• In order to achieve this, the authors discretize the simulation domain into m cells (they use the same cells as for the electrical conductivity problem) and to each cell they assign a source current density.
• For each cell (numbered from 1 to m), the authors 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 they compute the resulting distributions of the potential at each of the nφ voltage electrodes recording the secondary voltages.
• 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.

5.1 Experiment setup

• The experiment consists in mixing two materials, a clean sand and some burning coal.
• This water gives the saturated sand a resistivity of around 68 Ohm.m.
• These IP data were acquired like in a self-potential survey, that is, the authors inject current between two electrodes A and B and they measure the potential between all the remaining electrodes with respect to the reference electrode (electrode #1).
• 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.
• Then the secondary voltage decay is measured using ten time windows of different lengths (see Table 2).

5.2 Inversion

• The authors 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.
• As described in Section 2, the authors use this observed resistance data to recover the electrical conductivity spatial heterogeneities.
• 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.
• Recent advances and applications in complex resistivity measurements, Geophysics, 40(5), 851–864.

Figures

Figure 8. 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 7. 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 9. 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 16. 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 15. 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 3. 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 11. 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 2. 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 6. 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 14. 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 1. 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 12. 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. 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 10. 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.

Full text

HAL Id: hal-02324216
https://hal.archives-ouvertes.fr/hal-02324216
Submitted on 21 Oct 2019
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
3D time-domain induced polarization tomography: a
new approach based on a source current density
formulation
A. Soueid Ahmed, A. Revil
To cite this version:
A. Soueid Ahmed, A. Revil. 3D time-domain induced polarization tomography: a new approach based
on a source current density formulation. Geophysical Journal International, Oxford University Press
(OUP), 2017, 213, pp.244-260. �10.1093/gji/ggx547�. �hal-02324216�

Geophysical Journal International
Geophys. J. Int. (2018) 213, 244–260 doi: 10.1093/gji/ggx547
Advance Access publication 2017 December 18
GJI Marine geosciences and applied geophysics
3-D time-domain induced polarization tomography: a new approach
based on a source current density formulation
A. Soueid Ahmed
1,2
and A. Revil
1,2
1
Universit
´
e Grenoble Alpes, CNRS, IRD, IFSTTAR, ISTerre, F-38000 Grenoble, France. E-mail: andre.revil@univ-smb.fr
2
Universit
´
e Savoie Mont Blanc, ISTerre, F-73000 Chamb
´
ery, France
Accepted 2017 December 17. Received 2017 December 13; in original form 2017 April 17
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 polar-
ization (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 in-
verse 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 con-
ductivity 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 exper-
iment 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.
Key words: 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 polariza-
tion (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 salin-
ity, temperature, water content, and the cation exchange capacity
(CEC) of the material (e.g. Waxman & Smits 1968; Shainberg et al.
1980;Revilet 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 shut-
ting 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
244

The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society

Induced polarization tomography 245
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 explo-
ration 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 per-
formances (e.g. Zimmermann et al. 2008), IP has become a very
important method to investigate a broad spectrum of environmen-
tal applications. One can cite, for instance, the study of contami-
nants plumes (e.g. Olhoeft 1984, 1985, 1986; Vanhala et al. 1992,
Vanhala 1997; Slater & Lesmes 2002; Kemna et al. 2004;Wain-
wright et al. 2014), the characterization of permeability and pore
size distribution (e.g. Sturrock et al. 1999;Revilet al. 2015c;Joseph
et al. 2016;Ostermanet al. 2016), and recently coal seam fires de-
tection 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 re-
laxation 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 mod-
els such as the dynamic Stern layer model (Rosen & Saville 1991;
Rosen et al. 1993;Revil2012, 2013), these Cole–Cole parameters
can be interpreted in terms of textural and electrochemical proper-
ties 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 proper-
ties 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;Revilet al. 2015a,b;Mao&Revil2016;Maoet al.
2016;Revilet al. 2017a,b). This is due to the very strong IP re-
sponse 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
) inver ted 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 recov-
ered 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 res-
olution 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 accom-
plished 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 electri-
cal conductivity field by the chargeability. Following this approach,
an apparent chargeability is obtained by solving the Poisson’s equa-
tion 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 cor-
responding 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 demon-
strated that the (primary) current injection creates a secondary cur-
rent density J
S
(t) in the conductive ground. This secondary current
is related to the primary current J
p
through the chargeability evo-
lution 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 simi-
lar 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,fora
preliminary s tep 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 differ-
ence 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 in-
stance 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 in-
verse problem because retrieving the current density J
S
(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 prob-
lem does not require the use of an iterative process. This means that
notable computational time savings can be made for the tomogra-
phy of the intrinsic chargeability field. New instrumentations that
are massively multichannels such as the IRIS Full waver instrument

246 A. Soueid Ahmed and A. Revil
can operate the way we advocate: all the stations measure simultane-
ously 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 in-
jections 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. Con-
sidering the secondary voltages as pseudo self-potential data is also
correct from a physical point of view since these secondary volt-
ages 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 exper-
iment 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 in-
duction 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
∇·
(
σ
∞
∇ψ
)
+ I δ
3
(r − r
0
) = 0, (1)
where σ
∞
(in S m
−1
) denotes the 3-D high frequency electrical con-
ductivity field (i.e. the instantaneous conductivity of a medium
submitted to an electrical field, electromagnetic induction ef-
fects neglected, see Fig. 1), ψ (in V) is the electrical poten-
tial field generated by the injection of the current I (in A),
δ
3
(r − r
0
) = δ(x − x
0
) δ(y − y
0
) δ(z − z
0
) is the Dirac distribution,
x, y, z represent the space locations and x
0
, y
0
, z
0
are the spatial co-
ordinates of the current sources locations. Eq. (1) is subject to the
following boundary conditions
ψ = α on 
1
, (2)
σ
∞
∇ψ · n = β on 
2
, (3)
with 
1
∪ 
2
= ∂ where ∂ denotes the simulation domain
boundaries, n is the outward unit vector perpendicular to 
2
.In
case of α = 0andβ = 0, we refer to these boundary conditions as
homogenous Dirichlet and homogenous Neumann boundary condi-
tions, 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 mesh-
ing outside the area of interest and performing some benchmark
testing).
Once eq. (1) is solved, one can compute the apparent re-
sistivity associated with the instantaneous conductivity through
Figure 1. 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 T
on
(from –T to 0). The secondary voltage is measured
after the primary current is shut down (during T
off
). (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
.
ρ
a
= K ψ
MN
/I ,whereK 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 map-
ping operator F(.) associating the electrical potential ψ
σ
∞
(Fig. 1)
to the electrical conductivity σ
∞
,thatis,ψ
σ
∞
= F(σ
∞
). This po-
tential ψ
σ
∞
is instantaneously recorded when the current injection
is turned on. Similarly, a potential ψ
σ
0
= F(σ
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
M = 1 −
σ
0
σ
∞
=
σ
∞
− σ
0
σ
∞
, (4)
which is equivalent to
σ
0
= σ
∞
(
1 − M
)
. (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

Induced polarization tomography 247
discussed by Seigel (1949, 1959). The DC conductivity is obtained
with
ψ
σ
0
= F
(
σ
∞
(
1 − M
))
. (6)
Similarly, an apparent chargeability can be computed as
M
a
=
ψ
σ
∞
− ψ
σ
0
ψ
σ
∞
=
F
(
σ
∞
)
− F
(
σ
∞
(
1 − M
))
F
(
σ
∞
)
. (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
M
t
i
,t
i+1
=
1
ψ
0

t
i+1
t
i
ψ
MN
(
t
)
dt, (8)
where M
t
i
,t
i+1
is the partial chargeability measured during the time
window [t
i
, t
i+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:
M
t
i
,t
i+1
=
1
ψ
0
(
t
i+1
− t
i
)

t
i+1
t
i
ψ
MN
(
t
)
dt, (9)
where the primary voltage is given in V and the secondary voltage in
mV. Partial chargeability can be related to the apparent chargeability
M
a
through the approximation:
M
t
i
,t
i+1
≈ M
a
(
t
i+1
− t
i
)
. (10)
Note that this approximation is valid only under the condition that
t
i
, t
i+1
 τ (see Florsch et al. 2011). This equation can be used
to determine for each quadripole ABMN (A and B being the cur-
rent electrodes and M and N the potential electrodes), an apparent
chargeability M
a
. 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 J
p
= σ
0
E (at low frequen-
cies ∇×E = 0 and therefore E =−∇ψ) has been applied from
−T to time 0 (Fig. 1a). During the injection of the primary cur-
rent, 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),
J
S
(t) =−M(t )J
p
. (11)
Eq. (11) means that the dipole moment associated with the polar-
ization 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 com-
plex 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,
J = J
p
+ J
S
, (12)
∇·J = I δ
3
(r − r
0
). (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,
J =
[
1 − M(t)
]
J
p
, (14)
J =−
[
1 − M(t)
]
σ
∞
∇ψ. (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.
J =−σ
0
∇ψ, (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,
M(0) = 0, (17)
M(+∞) = M. (18)
For a Debye model, we have
M(t) = M

1 − exp

−
t
τ

, (19)
which is the classical behaviour for an RC circuit. For a Cole–Cole
model, we have
M(t) = MG

t
τ
; c

, (20)
G

t
τ
0
; c

= 1 −
∞

n=0
(−1)
n

t
τ
0

nc
(1 + nc)
, (21)
G

t
τ
0
; c

=
∞

n=1
(−1)
n+1

t
τ
0

nc
(1 + nc)
, (22)
where ( . ) denotes the Euler gamma function defined by,
(x) =

∞
0
u
x−1
e
−u
du, (23)
where x > 0. Note that the series development in eq. (22) converges
very slo wly for t > 10 τ and c < 1. It is easy to show that for c = 1,
we recover the Debye model,
G

t
τ
; c = 1

= 1 −
∞

n=0
(−1)
n

t
τ

n
(1 + n)
, (24)
G

t
τ
; c = 1

= 1 −
∞

n=0

−
t
τ

n
n!
, (25)

Frequently asked questions

Q1. What contributions have the authors mentioned in the paper "3d time-domain induced polarization tomography: a new approach based on a source current density formulation" ?

HAL this paper is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not.