# Spectral induced polarization porosimetry

A. Revil, Nicolas Florsch, Christian Camerlynck

► To cite this version:

A. Revil, Nicolas Florsch, Christian Camerlynck. Spectral induced polarization porosimetry. Geophysical Journal International, 2014, 198 (2), pp.1016-1033. 10.1093/gji/ggu180. hal-01196405

HAL Id: hal-01196405

https://hal.science/hal-01196405

Submitted on 9 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific 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 diffusion de documents
scientifiques 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spectral induced polarization porosimetry

A. Revil,<sup>1,2</sup> N. Florsch<sup>3,4</sup> and C. Camerlynck<sup>4</sup>

Accepted 2014 May 14. Received 2014 May 11; in original form 2014 January 24

## SUMMARY

Induced polarization is a geophysical method looking to image and interpret low-frequency polarization mechanisms occurring in porous media. Below 10 kHz, the quadrature conductivity of metal-free sandy and clayey materials exhibits a distribution of relaxation times, which can be related to the pore size distribution of these porous materials. When the polarization spectra are fitted with a Cole–Cole model, we first observe that the main relaxation time is controlled by the main pore size of the material and that the Cole–Cole exponent *c* is never much above 0.5, a value corresponding to a Warburg function. The complex conductivity is then obtained through a convolution product between the pore size distribution and such Warburg function. We also provide a way to recover the pore size distribution by performing a deconvolution of measured spectra using the Warburg function. A new dataset of mercury porosimetry and induced polarization data of six siliciclastic materials supports the hypothesis that the Cole–Cole relaxation time is strongly controlled by the pore size, and especially the characteristic pore size corresponding to the peak of the pore size distribution from mercury porosimetry. The distribution of the pore throat sizes of these materials seems fairly well recovered using the Warburg decomposition of the spectral induced polarization spectra but additional data will be needed to confirm this finding.

**Key words:** Electrical properties; Hydrogeophysics; Permeability and porosity.

## 1 INTRODUCTION

Induced polarization characterizes the ability of a porous material to store reversibly electrical charges at low frequencies (≤10 kHz). The geophysical method named from the measurement of these low-frequency polarization mechanism seems very promising in environmental geosciences and hydrogeology to characterize the permeability distribution of the subsurface (e.g. Kemna 2000; Binley et al. 2005; Hördt et al. 2007; Revil & Florsch 2010; Attwa & Günther 2013), to differentiate sedimentary formations (Attwa & Günther 2012; Gazoty et al. 2012), to monitor biodegradation processes (Davis et al. 2006; Revil et al. 2012c) and to delineate contaminant plumes and to monitor their remediation (e.g. Vanhala 1997; Deceuster et al. 2005; Schmutz et al. 2010; Vaudelet et al. 2011a,b; Chen et al. 2012; Deceuster & Kaufmann 2012; Schwartz & Furman 2012). Reviews of recent progress in the field of induced polarization applied to hydrogeophysical problems can be found in Kemna et al. (2012) and Revil et al. (2012a).

Low-frequency polarization is usually associated with the fact that the movement of ions in a porous material is governed not only by Coulombic forces but also by concentration gradients (both potentials can be combined to form electrochemical potentials, for example, Marshall & Madden 1959). Minerals are generally coated

by an electrical double layer formed by a Stern layer of sorbed counterions and a diffuse layer of ions attached to the mineral surface through Coulombic interaction with the net electrical charge on the surface of the minerals (Avena & de Pauli 1998). The existence of this electrical double layer and specially the Stern layer seems to have a strong influence on electrical polarization. This has been proven recently by changing the type of specific cations sorbed in the Stern layer and looking at the effect of these changes upon the magnitude of the polarization (see Vaudelet et al. 2011a,b).

Induced polarization measurements can be performed either in the frequency domain (spectral induced polarization, SIP) or in the time domain. In addition to a normalized chargeability (difference between the low and high-frequency electrical conductivities), a distribution of relaxation times can be inverted from induced polarization spectral data by performing a deconvolution of the response using a transfer function for the system (see Tong et al. 2006a,b; Tarasov & Titov 2007, for time-domain induced polarization and Kemna 2000; Ghorbani et al. 2007; Chen et al. 2008; Nordsiek & Weller 2008 for frequency-domain induced polarization). Using a mechanistic model, the distribution of relaxation times can be associated in turn to a distribution of polarization length scales, for instance a distribution of grain sizes or pore sizes (see Tong et al.

<footnote n="1"><sup>1</sup>Department of Geophysics, Colorado School of Mines, Golden, CO 80401, USA. E-mail: arevil@mines.edu</footnote>

<footnote n="2"><sup>2</sup>CNRS-UMR 5559-LGIT, Université de Savoie, Equipe volcan, F-73376 Le-Bourget-du-Lac, France</footnote>

<footnote n="3"><sup>3</sup>UMMISCO UMI n°209, 32 Av. Henry Varagnat, F-93143 Bondy cedex, France</footnote>

<footnote n="4"><sup>4</sup>Sorbonne Universités, UPMC Univ Paris 06, UMR7619 Métis, F-75005 Paris, France</footnote>
2004, 2006a,b; Revil & Florsch 2010; Revil et al. 2012b; Bücker & Hördt 2013).

In this paper, we first use a simple Cole–Cole model to describe the low-frequency complex conductivity of porous media. In this case, we will also provide evidences that the Cole–Cole exponent c is comprised between zero (for very broad pore size or particle size distributions) and 0.5 (in the case of very narrow pore size or particle size distributions). This observation points out that the correct transfer function that should be used to perform the deconvolution of the spectra is not the Debye function (a Cole–Cole model with a c-exponent of 1) but a Warburg function (a special case of the Cole–Cole function with a c-exponent equal to 0.5). Then, we will develop a method to perform the deconvolution of induced polarization spectra using the Warburg function and to analyse the distribution of relaxation times in terms of pore size distribution. This approach may offer a non-intrusive, non-destructive and environmentally friendly way to get porosimetry data from SIP data rather than using mercury porosimetry.

## 2 A COLE-COLE MODEL FOR THE COMPLEX CONDUCTIVITY

The goals of this section are (1) to demonstrate that the main relaxation time in a Cole–Cole model of complex conductivity is related to the main pore size of a porous material and (2) to show that for very narrow pore size distributions (close to a delta function), the Cole–Cole exponent is close to 0.5.

### 2.1 Definitions

In frequency-domain induced polarization, we record the magnitude of the electrical conductivity from the current and the measured voltage (corrected by a geometrical factor depending on the position of the electrodes and boundary conditions) and a phase between the current and the voltage in response to a periodic harmonic current. We note $\omega = 2 \pi f$ the angular frequency in rad s⁻¹, f is the frequency in Hz and $i = (-1)^{1/2}$ the pure imaginary number, AB the current electrodes and MN the voltage electrodes. The results can be expressed as a complex conductivity $\sigma^*$ or a complex resistivity $\rho^* = 1/\sigma^*$. The relationships between the modulus of this conductivity $|\sigma|$ and the phase $\varphi$ and the real and imaginary components of the conductivity, $\sigma'$ and $\sigma''$, are given by $\sigma^* = |\sigma| \exp(i\varphi) = \sigma' + i\sigma''$, with $|\sigma| = \sqrt{\sigma'^2 + \sigma''^2}$ and $\tan \varphi = \sigma''/\sigma'$. In Section 2.2, we provide a simple general model to describe the complex conductivity of clayey or clean sands and sandstones.

### 2.2 Model description

Fig. 1 depicts three polarization mechanisms associated with the polarization of a mineral grain coated by an electrical double layer coating the surface of the grains (e.g. Wang & Revil 2010). They include (1) the deformation of the diffuse layer, which occurs at very high frequency (>10 kHz, see the so-called Debye–Falkenhagen effect, Falkenhagen 1934), (2) the polarization of the Stern layer (Schwarz 1962; Leroy et al. 2008) and (3) the diffusion or membrane polarization, which is related to the difference in the transference number for the cations and anions during their migration in the porous material (Marshall & Madden 1959; Bücker & Hördt 2013). It seems that the main polarization mechanism is associated with the Stern layer since we still observe a strong polarization at

<figure><img src="image_2.png" /><figcaption>Figure 1. Polarization of a single grain. (a) Electrical double layer (Stern plus diffuse layers) at equilibrium. (b) Polarization mechanisms of the grain under an electrical field E. Three mechanisms occur: (i) the Falkenhagen contribution is due to the deformation of the diffuse layer and is characterized by a very fast relaxation time. (ii) The Stern layer polarization is due to the migration of the counterions in the Stern layer. (iii) The membrane polarization is due to the increase of salinity in the direction of the electrical field and the decrease of salinity on the other side of the grain (actually this increase/decrease of the salinity depends on the value of the transference numbers for the cations and anions).</figcaption></figure>

the isoelectric point of the mineral surface (zeta potential equals zero) for which there is no electrical diffuse layer (e.g. Revil et al. 2013a,b).

A very popular complex conductivity model is the Cole–Cole function (e.g. Tarasov & Titov 2013):

$$ \sigma^*(\omega) = \sigma_{\infty} - \frac{M_n}{1 + (i\omega\tau_0)^c}, \qquad (1) $$

where the normalized chargeability is traditionally defined by (e.g. Kemna 2000),

$$ M_n = \sigma_{\infty} - \sigma_0 \geq 0, \qquad (2) $$

while the chargeability is defined by $M = 1 - \sigma_0/\sigma_\infty$ (Kemna 2000), c denotes to the Cole–Cole exponent ($0 \le c \le 1$), $\tau_0$ denotes the characteristic relaxation time (or time constant) and $\sigma_0$ and $\sigma_\infty$ denote the low-frequency and high-frequency asymptotic limits of the electrical conductivity, respectively. This model was derived from the paper of Cole & Cole (1941) in which the Cole–Cole expression was derived to describe the complex permittivity of polar fluids. Florsch et al. (2012) noted that the Cole–Cole model in complex conductivity is not equivalent to the Cole–Cole model used by Pelton et al. (1978) to describe the complex resistivity of porous rocks but some relationships can be drawn between the two
models (see Florsch et al. 2012; Tarasov & Titov 2013, for further details). In the following, we will call ‘Warburg model’, a special case of the Cole–Cole model described above by eq. (1) with *c* fixed to 0.5. We will call ‘Debye model’, a special case of the Cole–Cole model described above with *c* fixed to 1.

In the frequency dependent conductivity model obtained through a volume-averaging approach by Revil (2013a,b), the low- and high-frequency conductivities entering eqs (1) and (2) are given by:

$$ \sigma_0 = \frac{1}{F} \sigma_w + \left(\frac{1}{F\phi}\right) \rho_S \beta_{(+)}(1-f) \text{CEC}, \quad (3) $$

$$ \sigma_{\infty} = \frac{1}{F} \sigma_w + \left(\frac{1}{F\phi}\right) \rho_S [\beta_{(+)}(1-f) + \beta_{(+)}^S f] \text{CEC}, \quad (4) $$

where *F* denotes the electrical formation factor (dimensionless), *φ* denotes the connected porosity (dimensionless), *σ*<sub>w</sub> (in S m<sup>-1</sup>) corresponds to the pore water conductivity, *f* (dimensionless) denotes the fraction of counterions in the Stern layer (typically ~0.90 for clays, see Revil 2013a,b), *ρ*<sub>S</sub> denotes the mass density of the solid phase (typically 2650 kg m<sup>-3</sup> for clay minerals and silica), *β*<sub>(+)</sub> corresponds to the mobility of the counterions in the diffuse layer [e.g. *β*<sub>(+)</sub>(Na<sup>+</sup>, 25 °C) = 5.2 × 10<sup>-8</sup> m<sup>2</sup>s<sup>-1</sup>V<sup>-1</sup>], *β*<sub>(+)</sub><sup>S</sup> denotes the mobility of the counterions in the Stern layer *β*<sub>(+)</sub><sup>S</sup>(Na<sup>+</sup>, 25 °C) = 1.5 × 10<sup>-10</sup> m<sup>2</sup>s<sup>-1</sup>V<sup>-1</sup> for clay minerals and *β*<sub>(+)</sub><sup>S</sup>(25 °C, Na<sup>+</sup>) = 5.2 × 10<sup>-8</sup> m<sup>2</sup>s<sup>-1</sup>V<sup>-1</sup> for silica grains, Revil 2013a,b) and CEC denotes the cation exchange capacity of the material (expressed in C kg<sup>-1</sup>). The cation exchange capacity denotes the total amount of cations that can get sorbed on the surface of a mineral at a given pH (generally ~7) and the product (1 − *f*)CEC denotes the amount of cations (per mass of grains) located in the diffuse layer whereas *f* CEC denotes the amount of cations (per mass of grains) located in the Stern layer.

## 2.3 Relation between the Λ-Parameter and the Cole–Cole relaxation time

Johnson *et al.* (1986) introduced a dynamic pore throat size Λ that can be determined from the distribution of the electrical potential in a pore network in absence of surface conductivity (e.g. Bernabé & Revil 1995). Revil *et al.* (2012b) found that the Cole–Cole relaxation time τ<sub>0</sub> entering eq. (1) seems to be controlled by this pore size Λ (rather than the grain size) according to,

$$ \tau_0 = \frac{\Lambda^2}{2D_{(+)}^S}, \quad (5) $$

where *D*<sub>(+)</sub><sup>S</sup> denotes the diffusion coefficient of the counterions in the Stern layer (the layer of sorbed counterions on the mineral surface, see Revil 2012, 2013a,b and Fig. 1).

Katz & Thompson (1987) developed a relationship between the permeability at saturation *k* and the percolation length scale *r*<sub>c</sub> using percolation principles *k* = *r*<sub>c</sub><sup>2</sup>/(226*F*). A comparison with *k* = Λ<sup>2</sup>/8*F* (Johnson *et al.* 1986; Avellaneda & Torquato 1991) yields *r*<sub>c</sub> ≈ 5.3 Λ. Using Λ<sup>2</sup> = 8*Fk* and Λ = *r*<sub>c</sub>/5.3 in eq. (5) yields:

$$ \tau_0 \approx \frac{r_c^2}{56 D_{(+)}^S}. \quad (6) $$

The diffusion coefficient entering into the expression of the Cole–Cole time constant (eqs 5 and 6) can be related to the mobility of the counterions in the Stern layer entering eqs (3) and (4), *β*<sub>(+)</sub><sup>S</sup>, by the Nernst–Einstein relationship *D*<sub>(+)</sub><sup>S</sup> = *k*<sub>b</sub>*T* *β*<sub>(+)</sub><sup>S</sup>/*q*<sub>(+)</sub>|,

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Phase spectrum for glass beads having the same size but different roughness. The low-frequency polarization is dominated by the response of the grains (suspensions) or the pores (dense porous material). The intermediate frequency range is clearly dominated by the effect associated with the roughness of the grains. The high frequency polarization is associated with dielectric effects (Maxwell–Wagner polarization including the high-frequency dielectric limit, see discussion in Revil 2013a). Data from Leroy *et al.* (2008), *C*<sub>f</sub> denotes the salinity of the pore water solution (NaCl) and *ω* the pulsation frequency.</figcaption></figure>

where |*q*<sub>(+)</sub>| is the charge of these counterions (e.g. |*q*<sub>(+)</sub>| = *e* for Na<sup>+</sup> where *e* denotes the elementary charge). The value *β*<sub>(+)</sub><sup>S</sup> (Na<sup>+</sup>, 25 °C) = 1.5 × 10<sup>-10</sup> m<sup>2</sup>s<sup>-1</sup>V<sup>-1</sup> yields *D*<sub>(+)</sub><sup>S</sup>(Na<sup>+</sup>, 25 °C) = 3.8 × 10<sup>-12</sup> m<sup>2</sup>s<sup>-1</sup> for clay minerals (see Revil 2013a,b). For the clean sands and sandstones, the mobility of the sodium in water [*β*<sub>(+)</sub>(Na<sup>+</sup>, 25 °C) = 5.2 × 10<sup>-8</sup> m<sup>2</sup>s<sup>-1</sup>V<sup>-1</sup>] leads to a diffusion coefficient of *D*<sub>(+)</sub><sup>S</sup>(Na<sup>+</sup>, 25 °C) = 1.32 × 10<sup>-9</sup> m<sup>2</sup>s<sup>-1</sup>.

## 2.4 Cole–Cole model versus experimental data

We compare now the previous model with experimental data. We first plot the phase lag for a pack of glass beads in the frequency range 3 mHz to 45 kHz (Fig. 2). These glass beads have identical particle sizes. The low-frequency polarization has previously associated either with grain sizes (Revil & Florsch 2010) or with pore sizes (Revil *et al.* 2012b; Bücker & Hördt 2013, this work). Note that the phase scales as *ω*<sup>1/2</sup> at low frequencies. We will show later that this behaviour is that predicted by a Warburg function. Fitting these data with a Cole–Cole model, we obtain a Cole–Cole exponent of 0.48 ± 0.3 (Fig. 2). The intermediate frequency range seems controlled partly by the polarization of heterogeneities associated with the roughness of grain surfaces (see Leroy *et al.* 2008). The phase at high frequencies scales as *ω* and this behaviour is a characteristic of the high frequency dielectric response −*i*ωε<sub>∞</sub> where ε<sub>∞</sub> denotes the high-frequency dielectric permittivity.

In Figs 3(a) and (b), we plot the complex conductivity spectrum of a fine sand characterized by a narrow grain size distribution with a mean of 100 µm. Despite the fact that the grain size distribution of this sand is narrow, it is not as narrow as for the glass beads used for the measurements displayed in Fig. 2. Fig. 3(b) shows that a portion of the spectrum is not well fitted by a Cole–Cole model because of the potential effect of the roughness of the grain-pore water interface between the main polarization frequency and the high-frequency dielectric effect (see Fig. 2). Fig. 3(c) shows a similar attempt to fit the data but without considering the data comprised between
<figure><img src="image_4.png" /><figcaption>**Figure 3.** In-phase and quadrature conductivities of a clean sand with a very narrow grain size distribution with a mean diameter of 100 µm (data from Vaudelet *et al.* 2011a). (a, b) For the quadrature conductivity, the data comprised between the low-frequency peak of the spectrum and the high-frequency dielectric effect is not fitted because of the potential effect of the roughness of the grain-pore water interface in this domain as shown in Fig. 2. $ε_0$ denotes the dielectric constant of vacuum. (c) Fit of the quadrature conductivities data from Vaudelet *et al.* (2011a) not taking into account the data potentially influenced by the roughness of the grains. Grain diameter 100 µm, formation factor 3.1 and permeability $2.85 \times 10^{-12} \text{ m}^2$. The pore size of the sand is roughly 8.4 µm. (d) Microphotograph of the sand.</figcaption></figure>

the peak frequency and the high frequency dielectric behaviour. The data are pretty well fitted with a Cole–Cole model with a $c$-exponent of $c = 0.39$, therefore much lower than the value $c = 1$ that would be associated with a Debye model and a bit smaller than the one found for the glass beads ($c = 0.48$). We also observe the high-frequency behaviour in $ω$ in Fig. 3 characterizing dielectric behaviour.

From the previous figures, it is possible that the broadness of the polarization is associated with the broadness of the grain or pore size distributions. We test now this idea using a broader database of experimental data. Fig. 4 shows the cementation exponent $m$ entering Archie's law $F = \phi^{-m}$ between the formation factor $F$ and the connected porosity $\phi$ (Archie 1942) and the Cole–Cole exponent $c$ distribution for a variety of sands and sandstones characterized by

different pore size distributions. Very sharp particle size distribution (PSD) corresponds to PSD distributed over less than one decade while very broad PSD corresponds to the PSDs spreading over three decades. For very well-sorted sands, $c$ is equal to $0.40 \pm 0.12$ (11 samples with the very sharp grain size distributions, see Fig. 4). This shows that the exponent $c$ in the Cole–Cole model does not take the full range of values between 0 and 1 but is restricted to the range 0 and 0.5. The upper bound, $c = 0.5$ is consistent with a Warburg model and not with a Debye model for which $c = 1$. When the distribution of pore sizes is very broad, the exponent $c$ has a tendency to be close to zero (flat distribution) while $m$ increases to values around 2. In conclusion, for very narrow grain size distribution (all the grains have essentially the same size), we would expect $c$ to be close to 1 according to the Debye model
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cross-plot of the cementation exponent ($m > 1$) entering Archie's law $F = \phi^{-m}$ in which $F$ is the intrinsic formation factor corrected for surface conductivity (Archie 1942) and the Cole–Cole exponent $c$ (according to the Warburg model discussed in the main text, $0 \le c \le 0.5$). There are two extreme end-regions denoted as A and B. Domain A is typical of well-sorted clean sands while domain B is typical of clay-rich materials. The plain line denotes the trend between $m$ and $c$. Data from Revil et al. (2013b, saprolites and clayey andstones) and Vinegar & Waxman (1984, samples 3477, 3336A, 3478, 101 and 102).</figcaption></figure>

but we actually find that $c$ is always smaller than ~0.50 as shown in Figs 3 and 4.

We discuss now the relationship between the Cole–Cole relaxation time and the main pore size for both clay-free and clayey granular porous materials using a broad database of published experimental data. We plot the Cole–Cole relaxation time versus the main pore size data for clean sands and sandstones (Fig. 5) and for shaly sands and sandstones (Fig. 6). We note that the trends are linear in a log–log scale with a power-law coefficient of two as predicted by eqs (5) and (6). The second observation is that the two trends are distinct indicating that the diffusion coefficient for the counterions in the Stern layer is not the same on the surface of pure silica and clays, in agreement with the observations made by Revil (2013a,b).

We fit the data shown in Fig. 5 with the Cole–Cole model (plain line) calculated using a mobility of $5.2 \times 10^{-8} \text{ m}^2\text{s}^{-1}\text{V}^{-1}$ (hence a diffusion coefficient $D_{(+)}^{S}(\text{Na}^{+}, 25^{\circ}\text{C}) = 1.32 \times 10^{-9} \text{ m}^2\text{s}^{-1}$ in eq. 6) for the counterions in the Stern layer equal to the mobility (diffusion coefficient) in water. There is a fair agreement between the model and the data. In Fig. 6, we perform the same comparison for clayey sands and sandstones. Again the model fits the data very well with the model computed with a diffusion coefficient of $D_{(+)}^{S}(\text{Na}^{+}, 25^{\circ}\text{C}) = 3.8 \times 10^{-12} \text{ m}^2\text{s}^{-1}$ determined from the ionic mobility inferred by Revil (2012, 2013a) for clay minerals (see Appendix A for the description of the in-phase and quadrature components of the complex conductivity). Our model indicates that the mobility of the counterions of clays is roughly 300 times smaller on the surface of clay minerals than the value obtained for silica sands (see discussion above at the end of Section 2.3).

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Cole–Cole relaxation time $\tau_0$ versus the pore size $\Lambda$ for clean sands. For the data from Koch *et al.* (2011), the pore size is determined from the median grain size and the formation factor using the relationship developed by Revil & Florsch (2010). Pore size is determined from mercury intrusion porosimetry (Binley *et al.* 2005), from $r_c = 5.6 (8 k/\phi^2)^{1/2}$ for the data of Tong *et al.* (2006a) and from permeability for the data of Vaudelet *et al.* (2011a,b). For the clean sands and sandstones, the diffusion coefficient is $D_{(+)}^{S}(\text{Na}^{+}, 25^{\circ}\text{C}) = 1.32 \times 10^{-9} \text{ m}^2\text{s}^{-1}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Cole–Cole relaxation time $\tau_0$ versus the pore size for clayey and shaly porous media. For the clayey material, the diffusion coefficient is estimated to be equal to $D_{(+)}^{S}(\text{Na}^{+}, 25^{\circ}\text{C}) = 3.8 \times 10^{-12} \text{ m}^2\text{s}^{-1}$ and is determined from the mobility of the counterions in the Stern layer (see Revil 2013 and discussion in the main text). Data from Scott & Barker (2003), Revil *et al.* (2013a), Comparon (2005), Lesmes & Frye (2001) and this work. The $\Lambda$-parameter is determined from the permeability and the formation factor using $\Lambda^2 = 8Fk$.</figcaption></figure>

# 3 ELECTRICAL POROSIMETRY OF SANDSTONES: THEORY

## 3.1 Rationale for a Warburg transfer function

Warburg (1899) developed the first solution to the diffusion equation with oscillating concentration as boundary condition. His approach
has been subsequently used in the literature to model a variety of polarization mechanisms. The complex conductivity of a Warburg model can be written as

$$ \sigma^* = \sigma_{\infty} - \frac{M_n}{1 + (i\omega\tau_0)^{1/2}}. \qquad (7) $$

Eq. (7) shows that the Warburg model is a special case of the Cole–Cole model with an exponent of $c = 0.5$. At low-frequencies, the phase or the quadrature conductivities scales as $(\omega)^{-1/2}$ as shown in Fig. 2 for glass beads with a uniform grain size. When decomposed with a Debye model, the (normalized) distribution of relaxation times for the Cole–Cole model is (Cole & Cole 1941)

$$ g(\tau) = \frac{1}{2\pi \tau} \frac{\sin[\pi(1-c)]}{\cosh[c \ln(\frac{\tau}{\tau_0})] - \cos[\pi(1-c)]}. \qquad (8) $$

In the special case of the Warburg model, we have therefore:

$$ g(\tau, c = 0.5) = \left(\frac{1}{2\pi\tau}\right) \frac{1}{\cosh\left[\frac{1}{2}\ln\left(\frac{\tau}{\tau_0}\right)\right]}, \qquad (9) $$

$$ g(\tau, c = 0.5) = \left(\frac{1}{\pi}\right) \frac{\left(\frac{\tau}{\tau_0}\right)^{1/2}}{\tau + \tau_0}. \qquad (10) $$

The idea we follow in the remaining part of the paper is that the Warburg function, rather than the Debye function, is the correct transfer function to determine the pore size or the pore size distribution. In the case that all the pores have the same size, a single pore size should be reflected by a single relaxation time, that is the Warburg function should be analysed in terms of a distribution of relaxation times $h(\tau)$ given by $h(\tau) = \delta(\tau - \tau_0)$ (i.e. a Dirac function).

The relationship between the distribution of relaxation times obtained with Debye and Warburg decompositions of a given complex conductivity spectrum is discussed further in Appendix B. Using the convolution theorem based on the superposition principle, the complex conductivity is given by,

$$ \sigma^* = \sigma_{\infty} - M_n \int_{0}^{\infty} \frac{h(\tau)}{1 + (i\omega\tau)^{1/2}} d\tau, \qquad (11) $$

The distribution of the intrinsic relaxation times, $h(\tau)$, corresponds to a normalized probability density and therefore obeys to

$$ \int_{0}^{\infty} h(\tau) d\tau = 1. \qquad (12) $$

We believe that there are two reasons to choose the Warburg function as transfer function rather than the Debye function:

(i) There is a group of theoretical transfer functions based on mechanistic models (Wong 1979; de Lima & Sharma 1992; Dukhin & Shilov 2002) describing the polarization of charged colloidal particles and granular materials (all with the same polarization length scale) with a transfer function close to a Warburg transfer function. The Warburg model is usually used to model a leaky capacitance. In our case, this would mean that the grains coated by the Stern layer behave as leaking capacitances. If we consider for instance a sand saturated by a NaCl solution, the controlling reaction for the interfacial impedance corresponds to the sorption/desorption of the sodium on the silica surface (Fig. 1b):

$$ > \text{SiO}^{-}\text{Na}^{+} \Leftrightarrow > \text{SiO}^{-} + \text{Na}^{+}, \qquad (13) $$

where ‘>’ refers to the crystalline framework and >SiO⁻ denotes the negative site at the surface of the silica SiO₂. It follows that

the polarization shown in Fig. 1 can be seen as a leaky capacitance characterized by a Warburg impedance model.

(ii) If we adopt a Debye decomposition of spectra to obtain relaxation times, the associated distribution of polarization length scale seems always much broader than the true distribution of polarization length scales (e.g. Leroy et al. 2008; Vaudelet et al. 2011a). This observation is consistent with the other observation made above that the spectral induced polarisation data of glass beads and silica sands fitted with a Cole–Cole model show a Cole–Cole exponent always smaller or equal to 0.50. In both cases, the Debye model is not the model describing the behaviour of a porous material with all the grains being the same and the observations point out that the Warburg model is the correct model to be used in deconvolving the complex conductivity data.

## 3.2 Relation between the relaxation time distribution (RTD) and the pore size distribution

We consider the case of a sandstone characterized by a (normalized) pore size probability distribution $f(r)$ where $r$ is a given pore radius. We assume that this pore size distribution $f(r)$ is associated with a (normalized) distribution of relaxation times $h(\tau)$. Assuming,

$$ \tau = \frac{r^2}{2D_{(+)}^S}, \qquad (14) $$

$(D_{(+)}^S)$ denotes the diffusion coefficient of the counterions), the two probability density functions $f(r)$ and $h(\tau)$ are related to each other by their probability distributions following the same type of analysis made by Revil & Florsch (2010). This yields,

$$ \tau h(\tau) = \frac{f(r)}{2F_h}, \qquad (15) $$

with,

$$ F_h = \int_{0}^{+\infty} f(r) d\ln r. \qquad (16) $$

The parameter $F_h$ can be seen as a normalization constant. Eqs (14) and (15) can be used to obtain an explicit function of the relaxation time using the distribution of pore sizes,

$$ h(\tau) = \left( \frac{D_{(+)}^S}{F_h} \right) \frac{f(r)}{r^2}. \qquad (17) $$

Eq. (17) can be used to compute the distribution of relaxation times from the distribution of pore sizes. Alternatively, pore size distribution can be constrained from the distribution of relaxation times inverted from complex conductivity spectra. In a second step, this distribution can be inverted to retrieve $f(r)$. Such an inversion would allow getting the pore size distribution from SIP data in the same way it is presently performed with nuclear magnetic resonance (NMR, see for instance Timur 1969; Weller et al. 2010). In the next section, we explain how we can perform a deconvolution of SIP spectra and use the result to determine the pore size distribution.

## 3.3 Warburg decomposition of SIP spectra

We consider below that a porous rock is a linear system in terms of the relationship between the pore size distribution and the resulting induced polarization spectra. We have seen in the previous sections that if the pore size distribution is described by a delta function, the associated RTD will be described by a delta function, and the
associated impulse function of the system is a Warburg function. Therefore the pore size distribution can be recovered by performing a deconvolution of the induced polarization spectra using this impulse or Green's function. We propose below a methodology to perform this task, which is based on a generalization of the approach proposed recently by Florsch et al. (2012) who used the Debye function to perform the deconvolution of complex conductivity spectra.

We start by using a complete expression of the effective complex conductivity derived above in eq. (11) and adding the dielectric polarization (electric displacement) to the complex conductivity:

$$ \sigma_{\text{eff}}^* = \sigma_{\infty} - M_n \int_{0}^{\infty} \frac{h(\tau)}{1 + (i\omega\tau)^{1/2}} d\tau + i\omega\epsilon_{\infty}, \quad (18) $$

In the following, the set {$M_n, h(\tau)$} is replaced by a non-normalized function $\tilde{h}(\tau)$ (to be inverted after discretization, see Florsch et al. 2012, Appendix A):

$$ \sigma_{\text{eff}}^* - i\omega\epsilon_{\infty} = \sigma_{\infty} - \int_{0}^{\infty} \frac{\tilde{h}(\tau)}{1 + (i\omega\tau)^{1/2}} d\tau. \quad (19) $$

Using the following transformation of variables (Florsch et al. 2012):

$$ \begin{cases} z = -\ln(\omega) \Leftrightarrow \omega = e^{-z} \\ s = \ln(\tau) \Leftrightarrow \tau = e^{s} \end{cases}, \quad (20) $$

Eq. (19) is written in a convolutive form:

$$ \sigma^* = \sigma_{\infty} - \int_{-\infty}^{\infty} G_s(s) \tilde{H}_z(z-s) ds. \quad (21) $$

with,

$$ \begin{cases} G_s(s) = \tau h(\tau) \\ \tilde{H}(z, s) = \frac{1}{1 + ie^{c(z-s)}} \end{cases}. \quad (22) $$

The discretization of eq. (21) on $N_G$ points leads to equations that corresponds to a Riemann sum with step $\Delta s$. By separating the real and imaginary parts, we obtain:

$$ \begin{cases} \sigma_k^{*\Re} = \sigma_\infty - \Delta s \sum_{j=1}^{N_G} \Phi_{kj}^{*\Re} G_j \\ \sigma_k^{*\Im} = -\Delta s \sum_{j=1}^{N_G} \Phi_{kj}^{*\Im} G_j + C e^{-z_k} \end{cases} \quad (23) $$

At this stage the system of equations is linear, but it should include also the fundamental property that the RTD is a probability function and must be necessarily positive. This property is insured by an additional change of variable:

$$ \forall j, \quad G_j = e^{G'_j}. \quad (24) $$

leading to the final model in $G'_j$:

$$ \begin{cases} \sigma_k^{*\Re} = \sigma_\infty - \Delta s \sum_{j=1}^{N_G} \Phi_{kj}^{*\Re} \cdot e^{G'_j} \\ \sigma_k^{*\Im} = -\Delta s \sum_{j=1}^{N_G} \Phi_{kj}^{*\Im} \cdot e^{G'_j} + C e^{-z_k} \end{cases}. \quad (25) $$

To solve that set of non-linear equations, we use the 'non-linear generalized least-square method' developed by Tarantola & Valette (1982). To optimize the damping factor devoted to smooth the solution of this ill-posed problem (see Tikhonov 1977), we use the L-curve method (see Florsch et al. 2012 for the application of the

L-curve to a very similar problem but based on the Debye model). We can also use a variant of the method proposed by Florsch et al. (2012) solving the problem in the Fourier space. A third possibility is to transform the result of the Debye decomposition into a Warburg decomposition distribution of relaxation times. This third idea is explored in Appendix B. In Appendix B, we derive a formula (eq. B17) that can be used to transform the distribution of relaxation time obtained from a Debye deconvolution to a distribution of relaxation time obtained from a Warburg deconvolution and vice versa. This formula can be very useful to researchers using existing codes based on the Debye deconvolution of the complex conductivity spectra and who want to transform their results in terms of a distribution of relaxation times using a Warburg deconvolution.

Once we have obtained the distribution of the relaxation times, we still need to determine the pore size distribution. We start with the relationship,

$$ \tau = \frac{r^2}{a D_{(+)}^S}, \quad (26) $$

where $r$ is a characteristic pore size, $a$ is a constant. It follows:

$$ r = \sqrt{\tau a D_{(+)}^S}, \quad (27) $$

Actually, $h(\tau)$ is not needed, only the logarithm of the distribution: $H(s) = \tau h(\tau)$ is required. The independent variable $s$ is given in $\log_{10}$ (we use the base 10 by choice). Therefore the transformation from $\tau$ to $r$ is just an affine transformation:

$$ \log_{10} r = \frac{1}{2} [\log_{10} \tau + \log_{10} (a D_{(+)}^S)]. \quad (28) $$

Therefore, we can use directly the RTD in $\log_{10} \tau$ to determine the distribution of pore sizes in $\log_{10} r$. This transformation corresponds only to a translation along the x-axis and a factor 1/2 applied to the scaling. A comparison between this procedure and experimental data is provided in Section 4.

# 4 ELECTRICAL POROSIMETRY: COMPARISON WITH EXPERIMENTAL DATA

## 4.1 Core samples

Six siliciclastic samples were collected from the Great Divide Basin (Wyoming, USA). Five samples are sandstones with angular to subrounded grains and sample #439 denotes a mudstone (Table 1). The samples have been chosen because they are characterized by a broad variability in pore sizes. Their permeability ranges from 3 to 1600 mD and their porosity ranges from 0.21 to 0.31 (Table 1). The formation factors were obtained through classical conductivity measurements at different salinities (see for instance Revil 2012) and were corrected for surface conductivity.

## 4.2 Complex conductivity measurements

The six samples were saturated under vacuum with the natural ground water saturating the formations from which they were extracted. The pore water conductivity of this ground water is 0.048 ± 0.013 S m⁻¹ at 25 °C and its composition is reported in Table 2. The complex conductivity measurements were taken using the ZEL-SIP04-V02 impedance meter at 25 frequencies from 1 mHz to 45 kHz (see Fig. 7 for the experimental setup). This impedance meter was built by Egon Zimmermann in Germany (Zimmermann
<table><caption>**Table 1.** Description of the core samples in terms of porosity $\phi$, permeability $k$, formation factor $F$, pore size $\Lambda = (8kF)^{1/2}$ and Cole–Cole relaxation time $\tau_0$ from the measured SIP spectra. The measured grain density is in the range 2596–2620 kg m$^{-3}$. All the samples are sandstones except Sample #439, which is a mudstone. 1 mD–10$^{-15}$ m$^2$. The pore sizes $r_c(50)$ and $r_c(\text{peak})$ denote the median of the pore size distribution and the value of the pore radius corresponding to the peak of the distribution, respectively.</caption><thead><tr><th>Sample</th><th>φ<br>(−)</th><th>k<br>(mD)</th><th>F<br>(−)</th><th>Λ<br>(μm) (1)</th><th>τ<sub>0</sub><br>(s) (2)</th><th>r<sub>c</sub>(50)<br>(μm) (3)</th><th>r<sub>c</sub>(peak)<br>(μm) (4)</th><th>τ<sub>0</sub><br>(s) (5)</th></tr></thead><tbody><tr><td>S499</td><td>0.265</td><td>1103</td><td>5.6 ± 0.4</td><td>7.0</td><td>2.6 ± 0.2</td><td>14.3</td><td>22</td><td>2.3</td></tr><tr><td>S498</td><td>0.206</td><td>35.9</td><td>9.0 ± 0.8</td><td>1.6</td><td>0.20 ± 0.01</td><td>1.98</td><td>4.5</td><td>0.10</td></tr><tr><td>S490</td><td>0.233</td><td>635</td><td>12.1 ± 0.4</td><td>7.8</td><td>3.2 ± 0.3</td><td>12.2</td><td>25</td><td>2.9</td></tr><tr><td>S493</td><td>0.232</td><td>115</td><td>18.3 ± 0.4</td><td>4.1</td><td>0.41 ± 0.04</td><td>4.4</td><td>6</td><td>0.17</td></tr><tr><td>S439</td><td>0.208</td><td>2.62</td><td>13.3 ± 0.7</td><td>0.53</td><td>0.023 ± 0.001</td><td>0.52</td><td>0.9</td><td>0.004</td></tr><tr><td>S436</td><td>0.306</td><td>1623</td><td>4.0 ± 0.3</td><td>7.2</td><td>25.5 ± 3.6</td><td>13.2</td><td>55</td><td>14.2</td></tr></tbody></table>

Notes: (1) Determined using $\Lambda^2 = 8Fk$ with $k$ determined from the capillary entry pressure (see Fig. 15); (2) using a Cole–Cole fit of the spectral induced polarization spectra; (3) median of the pore size distribution from Hg porosimetry; (4) peak of the pore size distribution and (5) predicted value using $\tau_0 = r_c(\text{peak})^2/(56D_{(+)}^S)$ with $D_{(+)}^S(\text{Na}^+, 25^{\circ}\text{C}) = 3.8 \times 10^{-12} \text{ m}^2\text{s}^{-1}$.

<table><caption>**Table 2.** Composition of the natural groundwater for the sandstones. TDS, total dissolved solids.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Value</th></tr></thead><tbody><tr><td>TDS</td><td>mg l<sup>-1</sup></td><td>318</td></tr><tr><td>Conductivity</td><td>µS cm<sup>-1</sup></td><td>479</td></tr><tr><td>pH</td><td>—</td><td>8.1</td></tr><tr><td>Alkalinity</td><td>mg l<sup>-1</sup></td><td>109</td></tr><tr><td>Na<sup>+</sup></td><td>mg l<sup>-1</sup></td><td>30.6</td></tr><tr><td>K<sup>+</sup></td><td>mg l<sup>-1</sup></td><td>3.9</td></tr><tr><td>Ca<sup>2+</sup></td><td>mg l<sup>-1</sup></td><td>65.0</td></tr><tr><td>Mg<sup>2+</sup></td><td>mg l<sup>-1</sup></td><td>3.1</td></tr><tr><td>Cl<sup>-</sup></td><td>mg l<sup>-1</sup></td><td>6.0</td></tr><tr><td>HCO<sub>3</sub><sup>-</sup></td><td>mg l<sup>-1</sup></td><td>123</td></tr><tr><td>SO<sub>4</sub><sup>2-</sup></td><td>mg l<sup>-1</sup></td><td>132</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Experimental setup. (a) Position of the current (A and B) and voltage (M and N) electrodes and sensitivity map (assuming the sample is homogeneous), which shows that there is a good sensitivity of the measurement over the cross-section area of the column. The sensitivity for the electrode array ABMN was computed in Vaudelet *et al.* (2011b). (b) ZEL-SIP04-V02 impedance meter build by Egon Zimmerman.</figcaption></figure>

*et al.* 2008a,b). The accuracy was ~0.1–0.3 mrad at frequencies <1 kHz. Ag–AgCl electrodes were used for both injection and potential electrodes and the electrode array was circumferential around the cylindrical core samples (Fig. 7). Benchmark tests of this equipment can be found in Revil & Skold (2011).

The spectra are shown in Fig. 8. The magnitude of the conductivity increases with frequency and the phase shows well-characterized

bell shapes. At the relaxation peak, the phase is in the range –6 (mudstone, #439) to –38 mrad (Sample #499). The peak of the phase lag occurs in the range 2 mHz (Sample #436) to 20 Hz (mudstone, #439). The electrical conductivities are in the range 0.02–0.12 S m$^{-1}$.

## 4.3 Pore size distribution

Mercury porosimetry measurements were performed on the six samples (Figs 9a–14a). Thin sections of the core samples are also shown in Figs 9(b)–14(b). The permeability reported in Table 1 was determined from the mercury intrusion experiments using a modified Swanson (1981) approach derived by Revil *et al.* (2014) based on the capillary entry pressure $p_e$ (in Pa),

$$p_e = \frac{2\gamma}{\sqrt{226}} \left( \frac{\phi}{\sqrt{k}} \right), \qquad (29)$$

where $\phi$ denotes the connected porosity, $\gamma$ represents the surface tension between water and air (71.99 ± 0.05) × 10$^{-3}$ N m$^{-1}$ and $k$ the permeability at saturation. We test the proportionality between the capillary entry pressure and the ratio $\phi/\sqrt{k}$ in Fig. 15 (the proportionally coefficient is fixed by the theory). Eq. (20) is used to predict the permeability from the porosity and the capillary entry pressure to Mercury for the different samples. When using the following formula $\Lambda^2 = 8Fk$ to determine the pore scale $\Lambda$ we obtain $\Lambda \approx (0.38\gamma)/p_e$, which is used to determine $\Lambda$ using the capillary entry pressure.

## 4.4 Data analysis

We found that the induced polarization spectra cannot be fitted by the Cole–Cole model. That said, the portion of the induced polarization data in the vicinity of the relaxation peak can be fitted very well by the Cole–Cole model described in Section 2.1 and Appendix A (using eqs A4 and A5). The corresponding fits are shown by Figs 9(c)–14(c) (the non-linear optimization is done with the least-square criterion using the Gauss–Newton method, see Chen *et al.* 2008). The values of the Cole–Cole parameters [normalized chargeability $M_n$, Cole–Cole exponent $c$ and Cole–Cole relaxation time $\tau_0$] are shown on the Figs 9(c)–14(c), see also Table 1 for the value of $\tau_0$].

In Fig. 16, we plot the Cole–Cole relaxation time with the pore scale $\Lambda$ (determined from the formation factor and the permeability,
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Phase and magnitude of the complex conductivity for the six samples investigated in this study. (a) Phase lag angle in mrad. (b) Magnitude of the conductivity in S m⁻¹. For the mudstone, the phase data at low and high frequencies are not shown because they were characterized by a high standard deviation. The lines are guides for the eyes.</figcaption></figure>

see Johnson *et al.* 1986), the median of the pore size distribution from mercury porosimetry and the value of the pore size corresponding to the peak of the pore size probability distribution. We see that the Cole–Cole relaxation time is best correlated (R² = 0.94) with the value of the peak of the pore size probability distribution. The relaxation time predicted by eq. (19) in Table 1 predicts fairly well the Cole–Cole relaxation time (R² = 0.89 in a linear–linear plot).

The results of the inversion of the spectra are reported in terms of relaxation times and pore size distributions (Figs 17–19). In these figures, we see both the results of the Debye and Warburg decomposition directly in terms of pore size distribution using the transform developed in Section 3.3. The Debye and Warburg decompositions fit the data equally well (Figs 17–19).

The pore size distributions resulting from the Debye decomposition (c = 1 in the equations of Section 3.3) are generally too broad with respect to the pore size distributions obtained from mercury

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Pore size distribution and quadrature conductivity. (a) Pore size distribution of Sample #499 (clayey sandstone). (b) Thin section. (c) The quadrature conductivity data (filled circles) was only fitted with a Cole–Cole model (plain line) for the frequencies around the peak frequency (R² = 0.99).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Pore size distribution and quadrature conductivity. (a) Pore size distribution of Sample 490 (clayey sandstone). (b) Thin section. (c) The quadrature conductivity data (filled circles) was only fitted with a Cole–Cole model (plain line) for the frequencies around the peak frequency (R² = 0.97).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Pore size distribution and quadrature conductivity. (a) Pore size distribution Sample #493 Clayey sandstone). (b) Thin section. (c) The quadrature conductivity data (filled circles) was only fitted with a Cole–Cole model (plain line) for the frequencies around the peak frequency (R² = 0.98).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 13.** Pore size distribution and quadrature conductivity. (a) Pore size distribution of Sample #436 (clayey sandstone). (b) Thin section. (c) The quadrature conductivity data (filled circles) was only fitted with a Cole–Cole model (plain line) for the frequencies around the peak frequency (R² = 0.93).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Pore size distribution and quadrature conductivity. (a) Pore size distribution of Sample #498 (clayey sandstone). (b) Thin section. (c) The quadrature conductivity data (filled circles) was only fitted with a Cole–Cole model (plain line) for the frequencies around the peak frequency (R² = 0.98).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Pore size distribution and quadrature conductivity. (a) Pore size distribution of Sample #439 (mudstone). (b) Thin section. (c) The quadrature conductivity data (filled circles) was only fitted with a Cole–Cole model (plain line) for the frequencies around the peak frequency (R² = 0.99).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Test of the relationship between the capillary entry pressure as a function of the porosity and permeability. The experimental data are from Huet *et al.* (2005). They correspond to 89 sets of mercury-injection (Hg-air) capillary pressure data. Core samples include both carbonate and sandstone lithologies. The permeability is here expressed in mD.</figcaption></figure>

porosimetry. The Warburg decomposition ($c = 0.5$ in the equations of Section 3.3) yields pore size distributions that are narrower and compatible with the observed pore size distribution from mercury porosimetry. The pore size distributions determined from Hg-porosimetry are shown in terms of the peak (with the vertical bar) and the width of the distribution at half the value of the peak (horizontal bar). The Warburg decomposition is generally reproducing the maximum and the width of the distribution from Hg porosimetry. That said, the details of the pore size distribution for small pore sizes are lost showing the limitations of the method.

The information on small pores seems lost due to the high-frequency dielectric behaviour.

Note that hat the Debye function is sharper in frequency domain and thus is expected to yield broader distributions in relaxation times and pore size and such broad pore size distributions in the pore sizes are not supported by the data. In Fig. 20, we compare the variance of the distributions determined from induced polarization porosimetry and mercury porosimetry. We see that SIP porosimetry is doing a fair job in reproducing the broadness of the pore size distributions, at least the main peak since the small pores cannot been observed through SIP porosimetry.

## 5 CONCLUSIONS

The following conclusions have been reached.

1. The distribution of relaxation times seems to be controlled by the distribution of the pore sizes. This observation is still not fully understood by the current mechanistic models. Pore scale numerical modelling may be needed to evaluate the ability of the local pore scale physics to reproduce these observations. The pore scale model developed recently by Bücker & Hördt (2013) may be a step in this direction.

2. For a collection of five sandstones and one mudstone core samples, we found that the Cole–Cole relaxation time is proportional to the characteristic pore size $\Lambda$ squared and to the pore size squared determined from the capillary entry pressure. The Cole–Cole relaxation time seems also correlated to the peak of the pore size distribution determined from Hg porosimetry.

3. The diffusion coefficients of the counterions moving along the Stern layer of pure silica sands and clays are consistent with the mobility of the Stern layer polarization (POLARIS) model developed by Revil (2012, 2013a,b). This shows the consistency in the model between the components describing the relaxation time and those needed to compute the amplitude of the quadrature conductivity at high and low frequencies. To the best of our knowledge, there is no other model able to reproduce such a consistency.

4. The variance as well as the peak of the pore throat size distribution can be recovered from the SIP data. Finer structures of

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Relationship between the Cole–Cole relaxation time and two different pore sizes reported in Table 1. (a) Correlation between the Cole–Cole relaxation time and the pore size $\Lambda$ determined using $\Lambda^2 = 8Fk$ ($R^2 = 0.90$) and $k$ determined from the capillary entry pressure. (b) Correlation between the Cole–Cole relaxation time and the peak of the pore size distribution determined from the mercury intrusion experiments ($R^2 = 0.94$). The best correlation coefficient is obtained with the pore size determined from the mercury porosimetry using the pore size corresponding to the peak of the distribution.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Fit of the SIP data of sample 498. (a) Fit of the amplitude (Debye decomposition). (b) Fit of the phase (Debye decomposition). (c) Distribution of pore throat sizes (Debye decomposition). (d) Fit of the amplitude (Warburg decomposition). (e) Fit of the phase with a Warburg decomposition. (f) Distribution of pore throat sizes (Warburg decomposition). (g, h, i, j, k, l) Same for Sample 493. The pore size distributions determined from Hg-porosimetry are shown in terms of the peak (with the vertical bar) and the width of the distribution at half the value of the peak (horizontal bar).</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Fit of the SIP data of sample 490. (a) Fit of the amplitude (Debye decomposition). (b) Fit of the phase (Debye decomposition). (c) Distribution of pore throat sizes (Debye decomposition). (d) Fit of the amplitude (Warburg decomposition). (e) Fit of the phase (Warburg decomposition). (f) Distribution of pore throat sizes (Warburg decomposition). (g, h, i, j, k, l) Same for Sample 439 (mudstone).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Sample 436</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Fit of the SIP data of sample 499. (a) Fit of the amplitude with a Debye decomposition. (b) Fit of the phase (Debye decomposition). (c) Distribution of pore throat sizes with the Debye decomposition. (d) Fit of the amplitude with the Warburg decomposition. (e) Fit of the phase (Warburg decomposition). (f) Distribution of pore throat sizes with the Warburg decomposition. (g, h, i, j, k, l) Same for Sample 436.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Comparison between the inverted pore size distribution using the SIP data (black filled circles) and the pore size distributions from the mercury porosimetry (peak and variance materialized by the vertical and horizontal bars, respectively). In these figures, we have moved the position of the peaks to match each other.</figcaption></figure>

the pore size distribution, especially for small pore sizes, seem more difficult to obtain from this procedure. The Debye decomposition seems unable to provide the correct variance regarding the pore size distribution for the six samples presented in this study. At the opposite, the Warburg decomposition provides narrower distributions of relaxation times, which translates to distributions of pore throat sizes that are compatible with those observed through Hg-porosimetry.

The present findings open the door to induced polarization porosimetry as currently done with NMR. NMR and SIP seem to be a natural tools to invert pore size distributions not only in the laboratory but also in the field through joint inversion based on a petrophysical model. Indeed both techniques are sensitive to the pore size and we have developed an approach to remove surface conductivity from the total measured conductivity of a porous rock by taking advantage of the relationship between surface conductivity and quadrature conductivity (Revil 2013b; Weller et al. 2013). In other words, this means we can isolate the formation factor and use the quadrature conductivity, corrected for the formation factor, to determine the pore size distribution. Therefore it is possible to formulate a joint inversion problem for the SIP and the NMR and take advantage of their very different sensitivity maps to get an image (tomogram) of the pore size distribution. If confirmed with additional datasets, this may open new doors in hydrogeophysics to image non-intrusively some characteristics of the pore size distribution with a combination of SIP/NMR data. Finally, drainage or immission of a porous material should have a specific signature on the distribution of relaxation times and this aspect will be investigated in a future contribution.

## ACKNOWLEDGEMENTS

We thank the National Science Foundation for the project 'SIP to invert permeability' (in collaboration with Lee Slater at Rutgers University and Andrew Binley at Lancaster University, NSF award NSF EAR-0711053) and the Office of Science (BER), US. Department of Energy (award DE-FG02-08ER646559). We thank Egon Zimmermann for the construction of his impedance meter, Markus Hilpert and two anonymous referees for their constructive comments regarding a previous version of our paper, Andreas Hördt, and an anonymous referee for their constructive comments regarding this version of the paper.

## REFERENCES

* Archie, G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics, *Trans. AIME*, **146**, 54–62.

* Attwa, M. & Günther, T., 2012. Application of spectral induced polarization (SIP) imaging for characterizing the near-surface geology: an environmental case study at Schillerslage, Germany, *Aust. J. Basic Appl. Sci.*, **6**(9), 693–701.

* Attwa, M. & Günther, T., 2013. Spectral induced polarization measurements for predicting the hydraulic conductivity in sandy aquifers, *Hydrol. Earth Syst. Sci.*, **17**, 4079–4094.

* Avellaneda, M. & Torquato, S., 1991. Rigorous link between fluid permeability, electrical conductivity, and relaxation times for transport in porous media, *Phys. Fluids A*, **3**, 2529–2540.

* Avena, M.J. & De Pauli, C.P., 1998. Proton adsorption and electrokinetics of an Argentinean Montmorillonite, *J. Colloid Interface Sci.*, **202**, 195–204.
Bernabé, Y. & Revil, A., 1995. Pore-scale heterogeneity, energy dissipation and the transport properties of rocks, *Geophys. Res. Lett.*, **22**(12), 1529–1552.

Binley, A., Slater, L.D., Fukes, M. & Cassiani, G., 2005. Relationship between spectral induced polarization and hydraulic properties of saturated and unsaturated sandstones, *Water Resour. Res.*, **41**, W12417, doi:10.1029/2005WR004202.

Börner, F.D., 1992. Complex conductivity measurements of reservoir properties, in *Proceedings of the Third European Core Analysis Symposium*, Paris, pp. 359–386.

Bücker, M. & Hördt, A., 2013. Analytical modelling of membrane polarization with explicit parametrization of pore radii and the electrical double layer, *Geophys. J. Int.*, **194**, 804–813.

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.

Chen, J., Hubbard, S.S., Williams, K.H., Flores Orozco, A. & Kemna, A., 2012. Estimating the spatiotemporal distribution of geochemical parameters associated with biostimulation using spectral induced polarization data and hierarchical Bayesian models, *Water Resour. Res.*, **48**, W05555, doi:10.1029/2011WR010992.

Cole, K.S. & Cole, R.H., 1941. Dispersion and absorption in dielectrics. I. Alternating current characteristics, *J. Chem. Phys.*, **9**, 341–351.

Componon, L., 2005. Etude expérimentale des propriétés électriques et diélectriques des matériels argileux consolidés, *PhD thesis*, Institut de Physique du Globe de Paris, 400 pp.

Davis, C., Atekwana, A., Atekwana, E., Slater, L.D., Rossbach, S. & Mormile, M.R., 2006. Microbial growth and biofilm formation in geologic media is detected with complex conductivity measurements, *Geophys. Res. Lett.*, **33**, L18403, doi:10.1029/2006GL027312.

Deceuster, J. & Kaufmann, O., 2012. Improving the delineation of hydrocarbon-impacted soils and water through induced polarization (IP) tomographies: A field study at an industrial waste land, *J. Contam. Hydrol.*, **136–137**, 25–42.

Deceuster, J., Chaballe, M. & Kaufmann, O., 2005. 3D resistivity and IP tomographies as efficient tools to monitor remedial actions over a gas plume, in *Proceedings of the Symposium on the Application of Geophysics to Engineering and Environmental Problems (SAGEEP)*, San Diego, CA, pp. 155–166.

de Lima, O.A.L. & Sharma, M.M., 1992, A generalized Maxwell–Wagner theory for membrane polarization in shaly sands, *Geophysics*, **57**, 431–440.

Dukhin, S.S. & Shilov, V.N., 2002. Nonequilibrium electric surface phenomena and extended electrokinetic characterization of particles, in *Electrokinetics and Electrophoresis, Surfactant Sci. Ser.*, Vol. 106, pp. 55–85, ed. Delgado, A.V., Marcel Dekker.

Fiandaca, G., Auken, E., Vest Christiansen, A. & 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.

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 Geophys.*, **10**, 517–531.

Gazoty, A., Fiandaca, G., Pedersen, J., Auken, E., Christiansen, A.V. & Pedersen, J. K., 2012. Application of time domain induced polarization to the mapping of lithotypes in a landfill site, *Hydrol. Earth Syst. Sci. Discuss.*, **9**, 983–1011.

Ghorbani, A., Camerlynck, C., Florsch, N., Cosenza, P., Tabbagh, A. & Revil, A., 2007. Bayesian inference of the Cole–Cole parameters from time and frequency domain induced polarization, *Geophys. Prospect*, **55**(4), 589–605.

Hördt, A., Blaschek, R., Kemna, A. & Zisser, N., 2007. Hydraulic conductivity estimation from induced polarization data at the field scale-the Krauthausen case history, *J. appl. Geophys.*, **62**, 33–46.

Huet, C.C., Rushing, J.A., Newsham, K.E. & Blasingame, T.A., 2005. Modified Purcell/Burdine model for estimating absolute permeability from

mercury-injection capillary pressure data, in *Proceedings if the International Petroleum Technology Conference*, Doha, Qatar, 21–23 November 2005, Paper IPTC 10994, 12 pp.

Johnson, D.L., Plona, T.J. & Kojima, H., 1986. Probing porous media with 1st sound, 2nd sound, 4th sound and 3rd sound, in *Physics and Chemistry of Porous Media*, Vol. 2, pp. 243–277, eds Jayanth, R., Banavar, J. & Winkler, K.W., Am. Inst. Phys. .

Katz, A.J. & Thompson, A.H., 1987. Prediction of rock electrical conductivity from mercury injection measurements, *J. geophys. Res.*, **92**(B1), 599–607.

Kemna, A., 2000. Tomographic inversion of complex resistivity, theory and application, *PhD thesis*, Bochum University, 196 pp.

Kemna, A. et al., 2012. An overview of the spectral induced polarization method for near-surface applications. *Near Surface Geophys.*, **10**, 453–468.

Koch, K., Kemna, A., Irving, J. & Holliger, K., 2011. Impact of changes in grain size and pore space on the hydraulic conductivity and spectral induced polarization response of sand, *Hydrol. Earth Syst. Sci.*, **15**, 1785–1794.

Leroy, P., Revil, A., Kemna, A., Cosenza, P. & Ghorbani, A., 2008. Spectral induced polarization of water-saturated packs of glass beads, *J. Colloid Interface Sci.*, **321**, 103–117.

Lesmes, D.P. & Frye, K.M., 2001. Influence of pore fluid chemistry on the complex conductivity and induced polarization responses of Berea sandstone, *J. geophys. Res.*, **106**(B3), 4079–4090.

Marshall, D.J. & Madden, T.R., 1959. Induced polarization, a study of its causes, *Geophysics*, **24**, 790–816.

Nordsiek, S. & Weller, A., 2008. A new approach to fitting induced polarization spectra, *Geophysics*, **73**(6), F235–F245.

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**, 588–609.

Revil, A., 2012. Spectral induced polarization of shaly sands: Influence of the electrical double layer, *Water Resour. Res.*, **48**, W02517, doi:10.1029/2011WR011260.

Revil, A., 2013a. Effective conductivity and permittivity of unsaturated porous materials in the frequency range 1 mHz–1GHz, *Water Resour. Res.*, **49**, 306–327.

Revil, A., 2013b. On charge accumulations in heterogeneous porous materials under the influence of an electrical field, *Geophysics*, **78**(4), D271–D291.

Revil, A. & Florsch, N., 2010. Determination of permeability from spectral induced polarization data in granular media, *Geophys. J. Int.*, **181**, 1480–1498.

Revil, A. & Skold, M., 2011. Salinity dependence of spectral induced polarization in sands and sandstones, *Geophys. J. Int.*, **187**, 813–824.

Revil, A., Karaoulis, M., Johnson, T. & Kemna, A., 2012a. Review: some low-frequency electrical methods for subsurface characterization and monitoring in hydrogeology, *Hydrogeol. J.*, **20**(4), 617–658.

Revil, A., Koch, K. & Holliger, K., 2012b. Is it the grain size or the characteristic pore size that controls the induced polarization relaxation time of clean sands and sandstones?, *Water Resour. Res.*, **48**, W05602, doi:10.1029/2011WR011561.

Revil, A., Atekwana, E., Zhang, C., Jardani, A. & Smith, S., 2012c. A new model for the spectral induced polarization signature of bacterial growth in porous media. *Water Resour. Res.*, **48**, W09545, doi:10.1029/2012WR011965.

Revil, A., Skold, M., Hubbard, S.S., Wu, Y., Watson, D. & Karaoulis, M., 2013a. Petrophysical properties of saprolites from the Oak Ridge Integrated Field Research Challenge site, Tennessee, *Geophysics*, **78**(1), D21–D40.

Revil, A., Eppehimer, J.D., Skold, M., Karaoulis, M., Godinez, L. & Prasad, M., 2013b. Low-frequency complex conductivity of sandy and clayey materials, *J. Coll. Interface Sci.*, **398**, 193–209.

Revil, A., Barnier, G., Karaoulis, M. & Sava, P., 2014. Seismoelectric coupling in unsaturated porous media: Theory, petrophysics, and saturation front localization using an electroacoustic approach, *Geophys. J. Int.*, **196**, 867–884.
Schmutz, M., Revil, A., Vaudelet, P., Batzle, M., Femenia Viñao, P. & Werkema, D.D., 2010. Influence of oil saturation upon spectral induced polarization of oil bearing sands, *Geophys. J. Int.*, **183**, 211–224.

Schwarz, G., 1962. A theory of the low-frequency dielectric dispersion of colloidal particles in electrolyte solution, *J. Phys. Chem.*, **66**, 2636–2642.

Schwartz, N. & Furman, A., 2012. Spectral induced polarization signature of soil contaminated by organic pollutant: Experiment and modeling, *J. geophys. Res.*, **117**, B10203, doi:10.1029/2012JB009543.

Scott, J. & Barker, R., 2003. Determining pore-throat size in Permo-Triassic sandstones from low-frequency electrical spectroscopy, *Geophys. Res. Lett.*, **30**(9), 1450, doi:10.1029/2003GL016951.

Swanson, B.F., 1981. A simple correlation between permeabilities and mercury capillary pressures, *J. Petrol. Tech.*, **33**(12), 2498–2504.

Tarantola, A. & Valette, B., 1982. Generalized nonlinear inverse problems solved using the least square criterion, *Rev. Geophys. Space Phys.*, **20**(2), 219–232.

Tarasov, A. & Titov, K., 2007. Relaxation time distribution from time domain induced polarization measurements, *Geophys. J. Int.*, **170**, 31–43.

Tarasov, A. & Titov, K., 2013. On the use of the Cole–Cole equations in spectral induced polarization, *Geophys. J. Int.*, **195**, 352–356.

Tikhonov, A.N., 1977. *Solutions of Ill Posed Problems*, Scripta Series in Mathematics, V.H. Winston.

Timur, A., 1969. Pulsed Nuclear Magnetic Resonance studies of porosity, movable fluid, and permeability of sandstones, *J. Petrol. Tech.*, **21**, 775–786.

Tong, M., Li, L., Wang, W. & Jiang, Y., 2006a. Determining capillary-pressure curve, pore size distribution and permeability from induced polarization of shaley sand, *Geophysics*, **71**, N33–N40.

Tong, M., Li, L., Wang, W. & Jiang, Y., 2006b. A time-domain induced-polarization method for estimating permeability in a shaly sand reservoir, *Geophys. Prospect.*, **54**, 623–631.

Tong, M., Wang, W., Li, L., Jiang, Y. & Shi, D., 2004. Estimation of permeability of shaly sand reservoir from induced polarization relaxation time spectra, *J. Petrol. Sci. Eng.*, **45**, 31–40.

Vanhala, H., 1997. Mapping oil-contaminated sand and till with the spectral induced polarization (SIP) method, *Geophys. Prospect.*, **45**, 303–326.

Vaudelet, P., Revil, A., Schmutz, M., Franceschi, M. & Bégassat, P., 2011a. Induced polarization signature of the presence of copper in saturated sands, *Water Resour. Res.*, **47**, W02526, doi:10.1029/2010WR009310.

Vaudelet, P., Revil, A., Schmutz, M., Franceschi, M. & Bégassat, P., 2011b. Changes in induced polarization associated with the sorption of sodium, lead, and zinc on silica sands, *J. Coll. Interface Sci.*, **360**, 739–752.

Vinegar, H.J. & Waxman, M.H., 1984. Induced polarization of shaly sands, *Geophysics*, **49**, 1267–1287.

Wang, M. & Revil, A., 2010. Electrochemical charge of silica surface at high ionic strength in narrow channels, *J. Coll. Interface Sci.*, **343**, 381–386.

Warburg, E., 1899. Über das Verhalten sogenannter unpolarisiebarer Elektroden gegen Wechselstrom (About the behaviour of so-called “impolarizable electrodes” in the present of alternating current), *Ann. Phys. Chem.*, **67**, 493–499.

Weller, A., Nordsiek, S. & Debschütz, W., 2010. Estimating permeability of sandstone samples by nuclear magnetic resonance and spectral-induced polarization, *Geophysics*, **75**(6), E215–E226.

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.

Wong, J., 1979. An electrochemical model of the induced-polarization phenomenon in disseminated sulfide ores, *Geophysics*, **44**(7), 1245–1265.

Zimmermann, E., Kemna, A., Berwix, J., Glaas, W., Munch, H.M & Huisman, J.A., 2008a. A high-accuracy impedance spectrometer for measuring sediments with low polarizability, *Measur. Sci. Technol.*, **19**(1–9), 105603.

Zimmermann, E., Kemna, A., Berwix, J., Glaas, W. & Vereeken, H., 2008b. EIT measurement system with high phase accuracy for the imaging of spectral induced polarization properties of soils and sediments, *Measur. Sci. Technol.*, **19**(1–9), 094010.

## APPENDIX A: THE COLE-COLE MODEL

The Cole–Cole model is usually used to fit complex conductivity spectra. For the complex conductivity, the real (in phase) and imaginary (quadrature) components of the complex conductivity are given by (Cole & Cole 1941),

$$ \sigma' = \sigma_{\infty} - \frac{1}{2} M_n \left\{ 1 - \frac{\sinh[c \ln(\omega \tau_0)]}{\cosh[c \ln(\omega \tau_0)] + \sin[\frac{\pi}{2}(1-c)]} \right\}, \quad (A1) $$

$$ \sigma'' = -\frac{1}{2} \frac{M_n \cos[\frac{\pi}{2}(1-c)]}{2 \cosh[c \ln(\omega \tau_0)] + \sin[\frac{\pi}{2}(1-c)]}, \quad (A2) $$

where $M_n = \sigma_{\infty} - \sigma_0 \ge 0$. Including the high frequency effect in the complex apparent conductivity (associated with the high-frequency displacement current) yields (Revil 2013a),

$$ \sigma_{\text{eff}}^* = \sigma_{\text{eff}}' - i\sigma_{\text{eff}}'' = \sigma^* - i\omega\varepsilon_{\infty}. \quad (A3) $$

The components of the apparent complex conductivity are,

$$ \sigma_{\text{eff}}' = \sigma_{\infty} - \frac{1}{2} M_n \left\{ 1 - \frac{\sinh[c \ln(\omega \tau_0)]}{\cosh[c \ln(\omega \tau_0)] + \sin[\frac{\pi}{2}(1-c)]} \right\}, \quad (A4) $$

$$ \sigma_{\text{eff}}'' = \frac{1}{2} \left\{ \frac{M_n \cos[\frac{\pi}{2}(1-c)]}{\cosh[c \ln(\omega \tau_0)] + \sin[\frac{\pi}{2}(1-c)]} \right\} + \omega \varepsilon_{\infty}. \quad (A5) $$

## APPENDIX B: CONNECTING THE DEBYE AND WARBURG DECOMPOSITIONS

The complex conductivity is written as:

$$ \sigma^*(\omega) = b + a \int_0^\infty \frac{h(\tau)}{1 + (i\omega\tau)^c} d\tau, \quad (B1) $$

where $h(\tau)$ is the function to be retrieved and $c = 0.5$ for a Warburg decomposition of the spectra. Making the classical change of variables, we obtain,

$$ H_s(s) = \tau h(\tau), \quad (B2) $$

$$ z = -\ln(\omega\tau_0) \Leftrightarrow \omega\tau_0 = e^{-z}, \quad (B3) $$

$$ s = \ln\left(\frac{\tau}{\tau_0}\right) \Leftrightarrow \tau = \tau_0 e^{\xi}, \quad (B4) $$

In eq. (B1), the function $h(\tau)$ is normalized. However, when performing the inversion, it is easier to retrieve a function encapsulating the constant $a$ and to normalize at the end. In the following, we can drop the constant $a$. A normalization of the retrieved $h(\tau)$ is required as soon as we want to retrieve this constant and a normalized distribution. In addition, we set, without loss of generality, $\tau_0 = 1$, and we get in the general case:

$$ \sigma_z(z) = b + \int_{-\infty}^{\infty} H_s(s) \Phi(z-s) ds, \quad (B5) $$

with,

$$ \Phi(z-s) = \frac{1}{1 + i^c e^{-c(z-s)}}. \quad (B6) $$

The convolution integral can be easily obtained through Fourier transform. Coming back to eq. (B1), which is for interest (here in
the case $c = 0.5$), we get $\omega^c = e^{-cz}$ and $\tau^c = e^{cs}$, in such a way we can re-write:

$$ \sigma_z(z) = b + \int_{-\infty}^{\infty} \frac{H_s(s)}{1 + i^c e^{-c(z-s)}} ds = \sigma'_z(z) + i\sigma''_z(z). \quad (B7) $$

Following here Florsch et al. (2012), we only consider the imaginary part in this Fourier approach. Calculating the imaginary part of $\sigma_z(z)$ is facilitated when using the notation $i^c = \exp(i\pi c/2)$. We get:

$$ \sigma_z''(z) = - \int_{-\infty}^{\infty} H_s(s) \frac{\sin\left(\frac{\pi}{2}c\right) e^{-c(z-s)}}{1 + 2\cos\left(\frac{\pi}{2}c\right) e^{-c(z-s)} + e^{-2c(z-s)}} ds. \quad (B8) $$

It follows, by using the Fourier-convolution theorem and defining in this step $\Omega_z(z)$:

$$ \begin{aligned} \text{FT}[\sigma_z''] &= -\text{FT}[H_s] \cdot \text{FT} \left[ \frac{\sin\left(\frac{\pi}{2}c\right) e^{-cz}}{1 + 2\cos\left(\frac{\pi}{2}c\right) e^{-cz} + e^{-2cz}} \right] \\ &= -\text{FT}[H_s] \cdot \text{FT}[\Omega_{z,c}]. \end{aligned} \quad (B9) $$

The Fourier transform in eq. (B8) can be determined numerically. We write eq. (B9) using the '∼' notation to denote the corresponding functions in the Fourier domain, $\eta$ is the variable in the Fourier dual space. We obtain $\tilde{\sigma}_{\eta}(\eta) = -\tilde{H}_{\eta}(\eta) \cdot \tilde{\Omega}_{\eta,c}(\eta)$, and therefore:

$$ H_s(s) = \text{FT}^{-1}[\tilde{H}_{\eta}(\eta)] = -\text{FT}^{-1}\left[\frac{\tilde{\sigma}_{\eta}(\eta)}{\tilde{\Omega}_{\eta,c}(\eta)}\right]. \quad (B10) $$

As explained in Florsch et al. (2012), this form is not expected to work well, because $\tilde{\Omega}_{\eta,c}(\eta)$ can be really close to zero and then having this function in the dominator will cause some instabilities in getting $H_s(s)$. Actually, one could assume that where $\tilde{\Omega}_{\eta,c}(\eta)$ is small, $\tilde{\sigma}_{\eta}(\eta)$ should be also small, making the ratio in eq. (B9) limited. It is not the case in practice, because the data are always noisy, and then $\tilde{\sigma}_{\eta}(\eta) = -\tilde{H}_{\eta}(\eta) \cdot \tilde{\Omega}_{\eta,c}(\eta)$ must be written as:

$$ \tilde{\sigma}_{\eta}''(\eta) = -\tilde{H}_{\eta}(\eta) \cdot \tilde{\Omega}_{\eta,c}(\eta) + N(\eta), \quad (B11) $$

where $N(\eta)$ denotes the non-modelled noise in the dual domain. Following Florsch et al. (2012), the Fourier inversion scheme can be stabilized by involving the Wiener regularization and a L-curve procedure to obtain the optimum damping. We write:

$$ H_s^{\text{Wiener}}(s) = \text{FT}^{-1}[\tilde{H}_{\eta}(\eta)] = -\text{FT}^{-1}\left[\frac{\tilde{\sigma}_{\eta}(\eta)\tilde{\Omega}_{\eta,c}^*(\eta)}{|\tilde{\Omega}_{\eta,c}(\eta)|^2 + \lambda^2}\right], \quad (B12) $$

where $\lambda$ is a damping factor, and (*) denotes for the complex conjugate of the function. The Warburg case corresponds to $c = 0.5$. Because the complex conductivity spectrum can be written as a convolution product, the Fourier transform appears to be an efficient tool to determine the kernel of the convolution integral. This also permits to swap between Warburg and Debye decompositions. We come back on eq. (B11), just adding the index $c$ to H to specify explicitly that we try to recover H for a chosen $c$. Then, the general equation is

$$ \tilde{\sigma}_{\eta}''(\eta) = -\tilde{H}_{\eta,c}(\eta) \cdot \tilde{\Omega}_{\eta,c}(\eta). \quad (B13) $$

The Debye and Warburg cases are given by:

$$ \tilde{\sigma}_{\eta}''(\eta) = -\tilde{H}_{\eta,1}(\eta) \cdot \tilde{\Omega}_{\eta,1}(\eta), \quad (B14) $$

$$ \tilde{\sigma}_{\eta}''(\eta) = -\tilde{H}_{\eta,\frac{1}{2}}(\eta) \cdot \tilde{\Omega}_{\eta,\frac{1}{2}}(\eta), \quad (B15) $$

respectively. Therefore, the Debye and Warburg cases are related to each other by:

$$ \tilde{H}_{\eta,1}(\eta) \cdot \tilde{\Omega}_{\eta,1}(\eta) = \tilde{H}_{\eta,\frac{1}{2}}(\eta) \cdot \tilde{\Omega}_{\eta,\frac{1}{2}}(\eta). \quad (B16) $$

Then if we have determined the RTD using a Debye decomposition, we can derive the corresponding RTD using a Warburg decomposition by computing:

$$ H_{z,\frac{1}{2}}(z) = \text{FT}^{-1}\left[\frac{\tilde{H}_{\eta,1}(\eta) \cdot \tilde{\Omega}_{\eta,1}(\eta)}{\tilde{\Omega}_{\eta,\frac{1}{2}}(\eta)}\right]. \quad (B17) $$

Conversely, from a Warburg to a Debye decomposition, we have:

$$ H_{z,1}(z) = \text{FT}^{-1}\left[\frac{\tilde{H}_{\eta,\frac{1}{2}}(\eta) \cdot \tilde{\Omega}_{\eta,\frac{1}{2}}(\eta)}{\tilde{\Omega}_{\eta,1}(\eta)}\right]. \quad (B18) $$

It is natural that these two transforms do not involve the data set, since it is just a transform between relaxation time distributions. Of course, this relationship can be generalized easily to the transform from any Cole–Cole RTD to another one if they only differ by their constant $c$.