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-
entic 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 diusion de documents
scientiques 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)
![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).](/figures/figure-8-relaxation-of-the-secondary-voltages-following-the-3keu72x5.webp)
![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.](/figures/figure-7-sketch-of-the-sandbox-the-tank-is-filled-with-water-2z2nv5oe.webp)
![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.](/figures/figure-3-geometry-of-the-synthetic-experiment-for-this-case-17zgyhpc.webp)