# Using the stretched exponential function for automatic processing of time-domain induced polarization data and further interpretation

Léa Lévy, H Che, A Weller

► To cite this version:

Léa Lévy, H Che, A Weller. Using the stretched exponential function for automatic processing of time-domain induced polarization data and further interpretation. Geophysical Journal International, 2026, 245 (1), pp.ggag010. 10.1093/gji/ggag010. insu-05532549

HAL Id: insu-05532549

https://insu.hal.science/insu-05532549v1

Submitted on 2 Mar 2026

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 CC BY 4.0 - Attribution - International License
# Using the stretched exponential function for automatic processing of time-domain induced polarization data and further interpretation

L. Lévy,<fnref n="1" />,<fnref n="2" /> H. Che and A. Weller <fnref n="3" />

Accepted 2025 December 26. Received 2025 December 2; in original form 2025 May 20

## SUMMARY

Time-domain induced polarization (TDIP) data carry spectral information that can be used for petrophysical interpretation. At the same time, TDIP data can be collected in the field more efficiently than frequency-domain induced polarization (FDIP) data, thanks to the use of square-wave signals. However, TDIP field data are prone to noise, particularly strong near industrial installations and urban areas, above conductive media and in cases where little current is injected. The integral chargeability is a useful parameter to smoothen out the signal but it precludes any spectral interpretation. Debye decomposition (DD) is recognized as one of the best methods for spectral interpretation but the extracted parameters are particularly affected by data noise. More generally, processing TDIP data before further analysis, such as inversion or spectral analysis, is usually necessary for any quantitative interpretation. We propose here an automatic processing algorithm, based on the Kohlrausch-Williams-Watts (KWW) function, which is very close to the Havriliak-Negami model in frequency-domain, that fulfills this need. The processing is completed by an empirical handling of early-time electromagnetic coupling effects to improve the overall performance. The resulting procedure, tested and validated on three data sets that cover a large range of contexts, electrode configurations and acquisition settings, is available as open-source MATLAB scripts. The proposed approach is especially useful for further extracting spectral information from TDIP data through DD. Thanks to the theoretical framework offered by the KWW function, the behaviour of the integral chargeability could be investigated in a systematic manner, using both synthetic and field TDIP data. Recommendations could be formulated on how to make use of the spectral information, while keeping the automatic processing transparent and accessible to unexperienced users. This work advances the use of TDIP in the field of environmental geophysics.

**Key words:** Electrical properties; Electromagnetic theory; Induced polarization; Downhole methods; Hydrogeophysics; Instrumental noise.

## 1 INTRODUCTION

Electrical resistivity and induced polarization methods have been used over more than a century (C. Schlumberger 1920) and their development is still ongoing (e.g. A. Kemna et al. 2012; P. Kessouri et al. 2019). Induced polarization (IP) measurements consist in applying an electric current, and measuring the transient voltage, caused by electric charges accumulation at interfaces between conducting and insulating media, in response to the current (C. Schlumberger 1920; J. Bertin & J. Loeb 1976; J.S.

Sumner 1976). The injected current can either have a square or a sinusoidal waveform. When square-wave currents are injected, the method is often called ‘time-domain’ IP (TDIP), while it is called ‘frequency-domain’ IP (FDIP) when a sinusoid is injected (e.g. A. Binley & L. Slater 2020). Recently, the potential of induced polarization methods for monitoring reactive processes has been illustrated by laboratory (e.g. A. Mellage et al. 2018; S. Izumoto et al. 2022; F. Rembert et al. 2023), and field studies (e.g. K.H. Williams et al. 2009; A. Flores Orozco et al. 2011; L. Lévy et al. 2024).

<footnote n="1">¹CNRS, UMR 6118 Géosciences Rennes, 35042, Rennes, France. E-mail: lea.levy@univ-rennes.fr</footnote>

<footnote n="2">²Engineering Geology, Lund Tekniska Högskola, Lund University, 221 000, Lund, Sweden</footnote>

<footnote n="3">³Institute of Geotechnology and Mineral Resources, Technische Universität Clausthal, 38678 Clausthal-Zellerfeld, Germany</footnote>
While most laboratory investigations of induced polarization rely on impedance spectra obtained with FDIP data, this approach is too time-consuming to be widely deployed in the field—even though a few cases have shown it is possible (e.g. A. Flores Orozco et al. 2011). Thanks to the use of square-wave signals instead of sinusoids, TDIP data can be collected in the field more efficiently than FDIP, while they also carry significant spectral information. Therefore, in field applications, TDIP data acquisition is usually preferred (e.g. T. Dahlin & V. Leroux 2012; G. Fiandaca et al. 2013). However, the interpretation of TDIP field data in terms of biogeochemical processes is hindered by several limits. In general, TDIP data are more prone to noise than resistivity data. Sources of noise include harmonic noise, self-potential drift and spikes, which can be handled to some extent by post-processing (P.I. Olsson et al. 2016). However, low-frequency erratic noise, particularly affecting the late times where the voltage signal level is at its lowest, cannot be handled by such post-processing (E. Zimmermann et al. 2008; T. Martin et al. 2020). In addition, electromagnetic (EM) coupling effects, which include inductive and capacitive coupling, have long been a major challenge in both FDIP and TDIP measurements. While advances in FDIP have been achieved through improved instrument design (T. Radic 2004; E. Zimmermann et al. 2008) and model-based corrections (H. Wang et al. 2021, 2024), approaches for mitigating EM noise in TDIP remain limited (A. Dey & H.F. Morrison 1973; P.K. Fullagar et al. 2000). In TDIP, strong EM noise, as well as low-pass digital filters affect early-time data, producing a characteristic kink in the decay curve when plotted in log-log space (J.S. He et al. 2008; G. Fiandaca et al. 2012; L.H. Cox et al. 2019).

Noise is usually stronger near industrial installations and urban areas. Noise is particularly problematic in conductive media (e.g. presence of clay minerals or saline pore water) where the maximum voltage resulting from the current injection can be quite low and therefore the voltage discharge after current turnover may rapidly enter the noise level. Noise is also problematic when little current is injected, either due to high contact resistance or to limitations of the transmitting unit. Noise in TDIP data can be handled by calculating the integral chargeability, which tends to smoothen out the voltage signal in a relatively effective way (L. Lévy et al. 2024). However, the integral chargeability poses three problems: (i) all spectral information are lost, (ii) the integral chargeability depends on the current waveform and in particular on the time-window for integration and (iii) there exists no petrophysical relation linking the TDIP integral chargeability to biogeochemical processes, since this parameter is very little studied in the laboratory. On the other hand, Debye decomposition (DD) is a relevant method for extracting spectral information from field TDIP data (A. Tarasov & K. Titov 2007), since there is an equivalence between the time- and frequency-domain DD parameters (T. Martin et al. 2021), providing certain conditions are met, such as a sufficient length of the current step. Therefore, the extracted DD parameters can be related to petrophysical properties studied in the laboratory such as the quantity of certain minerals (e.g. metallic particles) or the pore space topology. Nevertheless, the DD may act as a noise amplifier if the voltage signal is not properly cleaned (T. Martin et al. 2021; L. Lévy et al. 2024) as it allows overfitting the TDIP discharge curves. Therefore, in order to be able to efficiently extract the spectral information carried by TDIP data and convert it to meaningful petrophysical information, for example, through

DD, data processing is needed to smoothen out the noise, while keeping the spectral information contained in the signal.

TDIP data processing is usually done by manually checking the voltage decays and rejecting single gates that deviate from the expected shape, removing the distorted portion at early times considered to be caused by EM coupling effects, or rejecting a full decay in case the expected shape cannot be identified at all (A. Gazoty et al. 2013; P.K. Maurya et al. 2018; A.S. Barfod et al. 2021). The need for automatable filtering schemes has been recognized and an automatic processing procedure based on empirical power-law fitting has been proposed (A. Flores Orozco et al. 2018). Machine learning approaches have also been developed to learn from manual processing but they are deemed too site specific to be generally applicable: the neural network needs to be trained again at each site (A.S. Barfod et al. 2021). On the other hand, there exists an analytical formula, corresponding to the Kohlrausch-Williams-Watts (KWW) function (R. Kohlrausch 1854; G. Williams & D.C. Watts 1970), that describes the shape of TDIP decays expected by users during manual processing: the stretched exponential shape (F. Alvarez et al. 1991; S. Kang & D.W. Oldenburg 2019). Therefore, what a user manually does during TDIP data processing is mostly visually checking whether the decay indicates a stretched exponential (SE) behaviour. We thus propose here a new methodology that uses this function, embedded in a full processing procedure, to automatically discard outlier gates or full decays that do not correspond to this shape. In addition, an automatic detection method is proposed to identify the kink point characteristic of EM effect and remove the affected early-time gates.

## 2 THEORY AND BACKGROUND

### 2.1 The stretched exponential function

The stretched exponential (SE) function, also referred to as Kohlrausch-Williams-Watts (KWW) function (R. Kohlrausch 1854; G. Williams & D.C. Watts 1970), has a generic form that describes the voltage response to a current step, also called 'step response' (see eq. 1).

$$f(t) = \alpha e^{-(\frac{t}{\tau_{KWW}})^{\beta}}, \quad (1)$$

where $\alpha$, $\beta$ and $\tau_{KWW}$ are three empirical parameters. Since the KWW function is a superposition of uncoupled Debye processes (F. Alvarez et al. 1991), a large range of IP decays can be represented by the SE function. The fact that a single KWW function corresponds to several Debye elements is illustrated, for example, in L. Lévy et al. (2024, fig. 8). Therefore, the SE function has been widely used to describe the time-domain relaxation behaviour of glass-forming liquids, porphyry mineral deposits and complex systems (R. Kohlrausch 1854; G. Williams & D.C. Watts 1970; F. Alvarez et al. 1991; R. Hilfer 2002; S. Kang & D.W. Oldenburg 2019).

An empirical relation exists between the coefficients of the time-domain SE function and the frequency-domain Havriliak-Negami (HN) relaxation function (S. Havriliak & S. Negami 1966), as shown by, for example, F. Alvarez et al. (1991) and S. Holm (2020). The HN function, described by function g in eq. (2), is a generalization of the Pelton model with one more parameter (S. Havriliak & S. Negami 1966) and it corresponds to the 'Generalized Cole-Cole' model (K.S. Cole & R.H. Cole 1941;
<figure><img src="image_2.png" /><figcaption>Figure 1. Examples of SE functions for the same value of α = 500 mVV⁻¹ and varying β and τKWW. Upper panels correspond to three fixed τKWW values (50-500-2000 ms) with varying β displayed as different colours (from 0.1 to 0.9). Bottom panels correspond to three fixed β values (0.1-0.5-0.9) with varying τKWW displayed as different colours (from 50 to 5000 ms).</figcaption></figure>

D.W. Davidson & R.H. Cole 1951; W. Pelton et al. 1983).

$$g(\omega) = \frac{\alpha_{HN}}{(1 + (i\omega\tau_{HN})^{\lambda})^{\gamma}}, \quad (2)$$

where α_HN, τ_HN, λ and γ are four empirical parameters.

The relationship among the parameters of the SE and HN models is first suggested by the fact that both models yield an accurate description of experimental data acquired on the same synthetic samples (F. Alvarez et al. 1991). This relationship is not analytical, since the HN and the SE relaxation functions are not exactly Fourier transforms of each other. In particular, the relation between the time constants in both functions depends on the exponent β in the SE function. For example, for β close to 1, τ_HN ≃ τ_KWW, and for β = 0.1, τ_HN = 300τ_KWW. The SE model has one parameter less so it can always be transformed to a HN model, but the contrary does not hold. The product λγ is close to β. Varying decay shapes with different values for SE parameters β and τ_KWW are illustrated in Fig. 1.

The parameter α in eq. (1) is the origin f(t = 0) of the polarizability decay. It thus corresponds to the ratio between the secondary voltage V_S when current is turned off and the maximum primary voltage V DC during the current-on time. This ratio is the theoretical definition of the time-domain chargeability by J.S. Sumner (1976), impossible to measure in practice.

The SE model, given its close link with the HN model—which covers a wider range of relaxation behaviours than the Cole-Cole model—as well as its connection to the theoretical definition of the chargeability through α, has thus the potential to describe a large range of time-domain IP decays found in natural environments and, at the same time, to provide relevant parameters for quantitative analysis. However, the SE function needs to be convolved with the appropriate current waveform used in the field in order to extract intrinsic parameters of the investigated medium. Here, we use a simplified approach without such convolution, where the SE model is only used as a tool to exclude noisy data (single gates or full decays) based on a SE-fitting algorithm and improve subsequent spectral analysis.

## 2.2 Pelton model and Debye decomposition

Several models have been developed to describe how ions reversibly accumulate at the interface between pore fluid and metallic particles, causing polarization. One of the most widespread phenomenological relaxation models to interpret FDIP data is the Pelton model (W.H. Pelton et al. 1978; A. Tarasov & K. Titov 2013; A. Weller & L. Slater 2022), given by eq. (3).

$$\rho^*(f) = \rho_{DC} \left[ 1 - m \left( 1 - \frac{1}{1 + (i2\pi f \tau)^c} \right) \right] \quad (3)$$

with four fitting parameters: the DC resistivity ρ DC, the chargeability m, the Pelton time constant τ and the so-called Cole–Cole exponent c.

However, the Pelton model, presented in eq. (3), is mostly applied to spectra measured on samples with metallic particles (of infinite electrical conductivity) of a single grain size or a log-normal grain size distribution. Variations of this model for two grain size fractions exist but the fitting parameters are not directly comparable to each other, and even the two grain size fractions model may not represent the complexity of naturally occurring distributions of pyrite crystals (L. Lévy et al. 2019).

A more general approach to analysing FDIP data is the Debye decomposition (DD) (S. Nordsiek & A. Weller 2008). With this approach, the frequency-dependent complex electrical resistivity ρ*(f) describes an IP spectrum that can be fitted by a superposition of N Debye terms (Pelton models with c = 1), as shown in eq. (4).

$$\rho^*(f) = \rho_{DC} \left[ 1 - \sum_{k=1}^{N} m_k \left( 1 - \frac{1}{1 + i2\pi f \tau_k} \right) \right], \quad (4)$$

where m_k and τ_k are pairs of partial chargeability and relaxation time of a single Debye model. The function m_k = f(τ_k) is also referred to as relaxation time distribution (RTD). The N chargeability values m_k, which are related to a set of pre-defined relaxation times τ_k and the value of DC resistivity ρ DC are the
result of a fitting procedure described by S. Nordsiek & A. Weller (2008). Two parameters can be integrated from the DD to simplify the interpretations: the total chargeability $m_{tot}$ (eq. 5) and the mean relaxation time and $\tau_{mean}$ (eq. 6).

$$m_{tot} = \sum_{k=1}^{N} m_k \qquad (5)$$

$$\tau_{\text{mean}} = \exp \left( \frac{\sum_{k=1}^{N} m_k \ln(\tau_k)}{m_{\text{tot}}} \right). \qquad (6)$$

The DD approach has also been applied to TDIP data, yet with a convolution of the step response with the current waveform, leading to a nonlinear optimization procedure (A. Tarasov & K. Titov 2007). Since the DD allows more degrees of freedom than the SE fit, it is not considered suitable as a data processing tool, as it tends to fit noise and especially the early times affected by EM coupling or low-pass filters.

## 2.3 Integral chargeability

The most common way to convert polarizability decays into a parameter that can be inverted or simply represented as a function of depth is to integrate the decays and obtain the so-called 'integral chargeability', $M_{int}$ (W.M. Telford et al. 1990; D. Mao et al. 2016). It represents the area enclosed by the discharge curve, $V(t)$, and its zero asymptote, in a given time-window, $[t_1 : t_2]$, divided by the primary voltage $V_{DC}$ and the time-window width, as described by eq. (7) (J. Bertin & J. Loeb 1976; J.S. Sumner 1976).

$$M_{\text{int}} = \frac{1}{t_2 - t_1} \int_{t_1}^{t_2} \frac{V(t)}{V_{\text{DC}}} dt. \qquad (7)$$

Compared to existing approaches for analysing IP decays, such as DD or fitting a Cole–Cole model, the integral chargeability presents a series of advantages: (i) it is rapid to compute, (ii) it is a simple formula that anyone can implement and (iii) it averages out noise across the decays. However, questions arise when it comes to both conversion of $M_{int}$ to petrophysical parameters and comparison of $M_{int}$ across data sets using different current waveforms. This will be further investigated in this study, based on the framework provided by the SE function.

# 3 MATERIALS AND METHODS

## 3.1 Algorithm for automatic TDIP data processing

The general idea of the algorithm developed here is to fit with an SE function the normalized voltage signal, also called polarizability, to discriminate between IP signal and noise. The polarizability, noted $\eta(t)$ in T. Martin et al. (2021) and A. Tarasov & K. Titov (2007), is defined by the ratio between the decaying voltage during the off-time, $V_{decay}(t)$, and the maximum voltage reached during the on-time, $V_{DC}$:

$$\eta(t) = \frac{V_{\text{decay}}(t)}{V_{\text{DC}}} = \frac{1}{2} \left( \frac{V_{\text{decay, +}}(t)}{V_{\text{DC, +}}} + \frac{V_{\text{decay, -}}(t)}{V_{\text{DC, -}}} \right), \qquad (8)$$

where $V_{decay,+}(t)$ and $V_{decay,-}(t)$ are the voltage time-series during the off-time, following the positive and negative current polarities, in case of a 50 per cent duty cycle. These voltage time-series are obtained from the full waveform signal after performing

systematic signal processing steps, such as 50 Hz noise removal and re-gating, that allow optimizing the signal-to-noise ratio (see e.g. P.I. Olsson et al. 2015; L. Lévy et al. 2024). The data processing proposed here starts after all these steps have been carried out, using the polarizability time-series calculated from the re-gated voltage.

For each datapoint—one datapoint corresponding to a quadrupole of electrodes, the processing input is an array of the form $(t_i_{centre}, \eta_i)$, where $t_i_{centre}$ is the gate centre and $\eta_i$ the polarizability value at gate $i$, for $i = 1..n$, $n$ being the number of gates (typically around 30 gates per datapoint for modern instruments). In practice, a series of steps are necessary before and after the SE fit itself to apply the processing procedure.

Step #1a of the algorithm consists in removing negative polarizability values from the array.

Step #2 checks whether at least 5 gates are left in the array. This step is repeated after each step that results in potentially removing certain gates. In case the number of remaining gates drops below 5, the procedure excludes the whole IP time-series of the datapoint and moves on to the next datapoint. While the threshold can be adjusted, the goal of this step is to avoid keeping polarizability curves with too few gates

At step #3, a log transform is applied to the $(t, f(t))$ array, corresponding in practice to $(t_i_{centre}, \eta_i)_{i=1..n}$.

Steps #4 and #5 deal with the early-time kink. While this feature is often attributed to pronounced EM coupling, it can also result from the low-pass filter built in the instrument, which may induce a similar distortion at early times (G. Fiandaca et al. 2012; T. Dahlin & V. Leroux 2012; P.I. Olsson et al. 2015). At step #4, a polynomial fit is performed on the $(\ln(t), \ln[f(t)])$ array, using the MATLAB polyfit function and a degree of 5 which was found to give the best performance. A test on the norm of residuals decides whether (i) the fit is good enough or (ii) a new polyfit should be carried out after removing the main outlier gates. The test works with the threshold of one on the norm of residuals and a threshold of 0.2 on the relative residuals of each individual gate. When the final polyfit reaches a norm of residuals below one, the main outlier gates are anyway removed before proceeding to step #5. In practice, the thresholds chosen here combined with the high polynomial degree only exclude gates that are very far from the trend. Step #5 calculates the curvature of the polynomial fit and finds all the local extrema (both maxima and minima). Curvature has been used to determine the optimal regularization parameter by identifying the 'corner' of the L-curve in ill-posed inverse problems (Y. Li & D.W. Oldenburg 1999). The earliest local extremum is considered as the characteristic 'kink point' after which EM noise becomes less important. All gates before this local extremum are removed, provided that its amplitude is above a certain thresholds (in order to avoid removing gates when there is no kink point attributed to the EM noise at early time).

Step #6a fits the $(\ln(t), \ln[f(t)])$ array, containing the remaining gates after step #5, with the function $g(x)$, corresponding to the log transform of the SE function as described by eq. (9).

$$\begin{align*} f(t) &= ae^{-(t/b)^c} \\ \Leftrightarrow \ln[f(t)] &= \ln(a) - \left(\frac{t}{b}\right)^c = \ln(a) - e^{\ln(t) - \ln(b)} \\ \Leftrightarrow g(x) &= a_1 - e^{\alpha x - b_1}, \end{align*} \qquad (9)$$
where the correspondence between the different variables and parameters is described in eq. (10).

$$ 
\begin{aligned}
g(x) &= \ln[f(t)] \\
x &= \ln(t) \\
a_1 &= \ln(a) \Leftrightarrow a = e^{a_1} \\
b_1 &= c \cdot \ln(b) \Leftrightarrow b = e^{b_1/c}
\end{aligned} 
\quad (10) 
$$

The fitting procedure at step #6 uses MATLAB nonlinear least-squares fitting tool, with parameters $a_1$, $b_1$, $c$ constrained by the boundaries described in eq. (11).

$$ 
\begin{aligned}
0 &< a_1 < 10 \\
-\infty &< b_1 < +\infty \\
0.1 &< c_1 < 1.
\end{aligned} 
\quad (11) 
$$

The log transform allows stabilizing the fit (parameters varying in shorter intervals) and reducing its dependence to the starting point. The starting point is chosen heuristically, first with $a_1 = 1$, $b_1 = 1$, $c = 0.5$. If the fit results in $R^2 < 0.99$, a new starting point is chosen, and this is repeated twice.

Step #6b excludes gates considered as outliers based on their absolute residuals, calculated as the difference between the fit, $g(x_i)$, and the data, $\ln(\eta_i)$, corresponding in practice to $\|g(x_i) - \ln(\eta_i)\| < 0.2$. The goal here is twofold: (i) improve the fitting quality, in particular the $R^2$, so that the $R^2$ can be used as a discriminating variable to decide whether or not to keep the full polarizability curve in the data set and (ii) clean the array to improve subsequent spectral inversions. Step #6c repeats the SE fit after removing outliers.

Finally, step #7, which is optional, excludes quadrupoles where the latest SE fitted parameters are out of range. The acceptable ranges considered here are described in eq. (12).

$$ 
\begin{aligned}
\alpha &< 1000 \quad [\text{mV V}^{-1}] \\
0.2 &\le \beta < 1 \\
1 &< \tau_{\text{KWW}} \quad [\text{mV V}^{-1}].
\end{aligned} 
\quad (12) 
$$

All threshold values and acceptable ranges mentioned here have been adjusted following trial-error work on the different data sets presented in the next section and may still be improved. Future users may want to adjust these values (starting points, number of tests, $R^2$-threshold, residual thresholds etc.) to improve the overall outlier selection.

## 3.2 Evaluation of the processing outcomes on different data sets

The processing algorithm was tested on three independent sets of data, whose characteristics are summarized in Table 1. These data sets were acquired at different sites, with different instruments, electrode configurations and acquisition settings. More details about the sites, survey objectives and technical features can be found in L. Lévy et al. (2019, 2022, 2024), for the IP surveys at Nesjavellir (Iceland), Krafla (Iceland) and Hvedemarken (Denmark), respectively. Two out of three data sets (Krafla and Hvedemarken) were previously processed independently of this algorithm, which also allows a comparison of the outcomes.

The performance of the processing procedure is evaluated in several ways. First, for the three data sets, the percentage of datapoints left is counted at different steps of the procedure.

Then, for Krafla and Hvedemarken data sets, the new processing results are compared to the previously existing processing, by counting the number of datapoints that are accepted by both, rejected by both or accepted in one and rejected in the other. The same is also done for individual IP gates. Finally, for Krafla data set, 2-D DCIP spectral inversion results based on manual and SE-based processing are also compared, using the EEMverter software, with maximum phase angle parametrization (G. Fiandaca et al. 2024). In particular, the final data misfit of the inversion is compared, using the same misfit formula as presented in L. Lévy et al. (2019, eq. 4) and considering 2 per cent and 10 per cent relative standard deviation for DC data and IP data, respectively. In addition, the performance is qualitatively evaluated by manually checking the IP data to understand what type of curves or circumstances create a discrepancy with the previously existing processing.

## 3.3 Further analysis of the decays by Debye decomposition

Debye decomposition (DD) is carried out on processed data instead of raw data in order to overcome noise amplification (T. Martin et al. 2021). This effect is illustrated in Figs B1 and B2 of the appendix in L. Lévy et al. (2024).

Similarly to the method presented by L. Lévy et al. (2024), the polarizability time-series, $\eta(t_i)$, can either be recalculated from the SE function obtained by the processing algorithm and validated after step #7 or consists of the original data where only the remaining gates after step #7 are considered. The resulting polarizability time-series are fitted by MATLAB non-negative linear least-square tool, using the function presented in eq. (13) and the optimization problem expressed in eq. (14).

$$ \eta(t) = \sum_{k=1}^{N} m_k e^{-t/\tau_k} \quad (13) $$

$$ \min_{x \ge 0} \|Cx - d\|_2^2, \quad (14) $$

where $x$ is a vector with the $N$ $m_k$ model parameters, $d$ is a vector with the $n$ measured data $\eta(t_i)$ and $C$ is the linear multiplier defined as $C_{i,k} = e^{-t_i/\tau_k}$, where $t_i$ are the centre times of the $n$ gates $i = 1..n$ and $\tau_k$ are the pre-defined $N$ relaxation times $k = 1..N$. During the procedure, most of the coefficients $m_k$ are kept to 0, due to the non-negative constraints. For Nesjavellir data set, the $\tau_k(k = 1..N)$ are defined as a set of $N = 300$ logarithmically spaced values of relaxation times in the range 3 – 2000 ms, covering approximately three decades of decay time. The fit eventually results in 10–20 non-zero coefficients.

Examples of Debye decomposition, shown in Fig. 2, illustrate that applying the DD directly to raw data either tends to overfit noise at early times (Figs 2a and c) or fails (Figs 2e and g). That's why we consider here that the SE-based processing is an important step before trying to extract spectral information from the IP decays. Moreover, it can be observed on Figs 2(b), (d) and (f), (h) that the SE function corresponds to a superposition of several Debye elements, indicating that the spectral information is not lost after the transformation.

It is important to note that the DD carried out here provides only qualitative parameters, albeit containing spectral information at a limited computation cost, because the current waveform is not taken into account in the procedure. More complex procedures involving the convolution of the Debye step response
<table><caption>**Table 1.** Description of the three data sets used to test the performance of the processing algorithm based on SE fitting.</caption><thead><tr><th></th><th>Nesjavellir</th><th>Krafla</th><th>Hvedemarken</th></tr></thead><tbody><tr><td>TDIP instrument</td><td>QL40-IP by ALT</td><td>Terrameter LS1 by ABEM</td><td>Terrameter LS2 by ABEM</td></tr><tr><td>Type of layout</td><td>Wireline logging</td><td>Surface DCIP</td><td>Cross-borehole DCIP</td></tr><tr><td>Electrode configuration</td><td>Pole-pole</td><td>Multi-gradient (current and voltage electrodes on separate, parallel, cables)</td><td>Single-borehole (dipole-dipole and multi-gradient), equatorial and parallel cross-borehole</td></tr><tr><td>Electrode spacing</td><td>0.4 and 1.6 m (16 and 64 inches)</td><td>40 m along a cable, 1 m between cables</td><td>2.5 m between boreholes, 30 cm along a borehole</td></tr><tr><td>Spatial resolution</td><td>25 cm vertically</td><td>40 m laterally, decreasing resolution with depth</td><td>0.15 × 0.5 m cells (vertical x lateral)</td></tr><tr><td>Spatial coverage</td><td>400 m vertically, one-dimensional</td><td>1.2 km laterally</td><td>10 m vertically, 8 m laterally</td></tr><tr><td>Type of electrodes</td><td>Embedded electrodes on probe</td><td>20-cm long sticks</td><td>Ring electrodes around PVC tube</td></tr><tr><td>Material for electrodes</td><td>Stainless steel</td><td>Stainless steel</td><td>Stainless steel</td></tr><tr><td>Electrode contact with subsurface</td><td>Probe moving in water-filled open hole</td><td>Electrodes digged in top soil, enhanced contact with saltwater and bentonite</td><td>Unconsolidated sediments collapsed around electrodes, saturated aquifer</td></tr><tr><td>Current waveform</td><td>50 per cent duty cycle with no stacking and 2-s on-time and off-time</td><td>50 per cent duty cycle with variable stacking as well as on and off times: on-time in the range 0.4–4 s; off-time in the range 1.3–8 s.</td><td>100 per cent duty cycle with no stacking and 2-s on-time (and 0-s off-time)</td></tr><tr><td>Voltage sampling rate (Hz)</td><td>450</td><td>1000</td><td>3750</td></tr><tr><td>Number of sub-data sets</td><td>Eight (Two wells measured four times each in the context of monitoring)</td><td>Seven profiles measured once</td><td>Three (same 4-borehole profile measured three times in the context of monitoring)</td></tr><tr><td>Total number of datapoints</td><td>28636</td><td>8315</td><td>6031</td></tr></tbody></table>

(eq. 13) with the actual pulse response used in the field would be necessary to quantitatively interpret the Debye parameters (A. Tarasov & K. Titov 2007; G. Fiandaca et al. 2012, 2013). Nevertheless, the same current waveform is used for all the Nesjavellir data, allowing meaningful comparison of the spectral parameters obtained. In addition, a comparison is carried out between SE and DD parameters, regardless of the current waveform.

# 4 RESULTS

## 4.1 Processing outcomes on the three data sets with the automatic outlier removal

This section illustrates how the outlier removal procedure works (Fig. 3) and summarizes the outcomes in terms of data kept versus data discarded for the three data sets (Table 2): Krafla, Hvedemarken and Nesjavellir.

Fig. 3 shows the different steps of the processing algorithm: (i) separating the positive (blue) and negative (red) gates and excluding the negative ones before proceeding (which could be reconsidered in future work), (ii) fitting the positive gates with a fifth-degree polynomial (blue line), (iii) calculating the curvature of the polyfit (red line) and identifying local extrema (numbers on the red line), (iv) carrying out a first SE fit on the remaining gates (green dots), after excluding outliers based on the red curvature, (v) if needed, carrying out a second SE fit (black line) on the remaining gates (black circles), after excluding outliers based on the first SE fit. The left column (Figs 3a, d, g) shows how the algorithm performs on IP curves with relatively high signal-to-noise ratio (SNR). The middle column illustrates cases of strong noise at early times (Figs 3b,e, most likely EM noise since the low-pass filter used by the Terrameter LS is at 1500 Hz) or erratic noise possibly related to a poor handling of the 50 Hz noise (Fig. 3h). The right column shows quadrupoles with low

SNR that are excluded from the data set by the automatic processing, because (i) too few gates are left to carry out the first SE fit (Fig. 3c), (ii) the first SE fit does not meet the $R^2 < 0.99$ criterion and too few gates are left after excluding outliers to carry out the second SE fit (Fig. 3f) and (iii) the signal cannot be fitted by a SE function due to high levels of noise (Fig. 3i).

Table 2 summarizes the removal rates of the automatic procedure. After completing step #6 of the processing, 78 per cent, 51 per cent and 66 per cent of the quadrupoles are left in the Nesjavellir, Krafla and Hvedemarken data sets, respectively. All these quadrupoles are fitted by a SE function with $R^2 > 0.99$. The additional step #7 checks whether the SE parameters belong to a pre-determined range, considered as reasonable based on a synthetic analyses not presented here (eq. 12). Typically, very flat SE models with $\beta < 0.2$ are too uncertain to fit given the limited time-range (Fig. 1). However, the ranges proposed in (eq. 12) can be discussed and future work may help better determine the range of parameters that can be fitted by spectral IP inversion softwares.

Next, the distribution of fitted $\alpha$, $\beta$ and $\tau_{KWW}$ values is analysed. Similar trends are found for the three data sets (Fig. 4) despite the variable current waveforms. Nevertheless, some variability is observed and further analysed here. Since Nesjavellir data set contains a number of datapoints six times larger than Krafla and Hvedemarken data sets after step #6 (Table 2), the model parameters fitted for this data set are sub-divided into separate plots for 16-inch and 64-inch electrode spacing and into different colours for the three wells monitored at the site (NN3, NN4 and NN7). The larger electrode spacing for the 64-inch configuration results in a smaller number of datapoints for each well, compared to the 16-inch configuration. In addition, NN4 is a shallower well (open hole in the depth-range 200–400 m) than NN3 and NN7 (open holes in the depth-range 200–550 m) and has thus a smaller number of datapoints.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Examples of Debye decomposition (DD) on raw data (left column) and SE fitted function (right column) for one quadrupole in Nesjavellir data set (a–d) and another in Krafla data set (e–h). The filled circles correspond to the raw data in (a) and (e) and to the SE fit in (b) and (f). The black line in (a)–(b) and (e)–(f) corresponds to the DD fit.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Examples of processed data illustrating the behaviour for different types of signals: clean IP signals, EM effect at early times or erratic noise. (a) Krafla, profile ISL-10, datapoint #1656; (b) Krafla, profile ISL-11, datapoint #1413; (c) Krafla, profile ISL-1, datapoint #349; (d) Hvedemarken, round 4, datapoint #304; (e) Hvedemarken, round 3, datapoint #1005; (f) Hvedemarken, round 5, datapoint #1725; (g) Nesjavellir, NN3, 64-inch spacing, Baseline, Depth = 328.53 m; (h) Nesjavellir, NN3, 64-inch spacing, +40 d, Depth = 203.40 m; (i) Nesjavellir, NN3, 64-inch spacing, +40 d, Depth = 380.65 m. For each column, the top, middle and bottom row correspond to selected decay examples of Krafla, Hvedemarken and Nesjavellir data sets, respectively. The blue and red dots correspond to the raw data, either positive or absolute negative, respectively. The blue hased line corresponds to the polyfit. The red line is the curvature of the polyfit. The black line corresponds to the SE fit. The green and black circles correspond to the remaining gates kept for the first and second SE fit, respectively.</figcaption></figure>

<table><caption>**Table 2.** Summary of fitting results with the stretched exponential function for the three data sets. The number of datapoints is indicated by '#dp'.</caption><thead><tr><th></th><th></th><th>Nesjavellir</th><th>Krafla</th><th>Hvedemarken</th></tr></thead><tbody><tr><td>Total</td><td>#dp</td><td>29418</td><td>8315</td><td>6031</td></tr><tr><td>Before step#6</td><td>#dp</td><td>24187</td><td>6068</td><td>4818</td></tr><tr><td></td><td>per cent</td><td>82 per cent</td><td>73 per cent</td><td>80 per cent</td></tr><tr><td>After step #6</td><td>#dp</td><td>22948</td><td>4221</td><td>4004</td></tr><tr><td></td><td>per cent tot</td><td>78 per cent</td><td>51 per cent</td><td>66 per cent</td></tr><tr><td>After step #7</td><td>#dp</td><td>21945</td><td>3198</td><td>3854</td></tr><tr><td></td><td>per cent tot</td><td>75 per cent</td><td>38 per cent</td><td>64 per cent</td></tr></tbody></table>

At Nesjavellir, SE models corresponding to 16-inch electrode spacing data show a double Gaussian distribution for $\alpha$ (Fig. 4a). The main peak, at $\alpha_{mean} = 50 mV V^{-1}$, is similar to the only peak observed for the Krafla (Fig. 4c). This peak is only observed for NN3 and NN7 wells, not for NN4, and is absent from the SE models corresponding to 64-inch data, for all wells (Fig. 4b). The second peak observed in Fig. 4(a), at $\alpha_{mean} = 150 mV V^{-1}$, is similar to the second peak observed for the SE models corresponding to 64-inch data (Fig. 4b) in all wells. While it is beyond the scope of this study to interpret these observations, it is interesting to note that most values found in this natural environment fall in the range 20–200 mV V$^{-1}$.

The distribution of $\beta$ values is much less variable across data sets. All wells show a peak at $\beta_{mean} = 0.3$ for Nesjavellir data, both for 16-inch (Fig. 4e) and 64-inch electrode spacing (Fig. 4f),
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Distribution of $\alpha$ (a–d), $\beta$ (e–h) and $\tau_{KWW}$ (i–l) values for SE models fitted to the different data sets. The four rows correspond to different data sets: Nesjavellir with 16-inch electrode spacing; Nesjavellir with 64-inch electrode spacing; Krafla; Hvedemarken. The different colours on Nesjavellir plots correspond to the three wells monitored at Nesjavellir. SE models with $\beta < 0.2$ or $\tau_{KWW} < 1$ ms are considered too uncertain.</figcaption></figure>

<table><caption>**Table 3.** Comparison of manual and automatic processing for Krafla and Hvedemarken data sets, in terms of quadrupoles as well as IP gates used in the inversion. Note that in Krafla data set, one electrode quadrupole usually corresponds to several datapoints because different pulse sequences were used for each electrode configuration.</caption><thead><tr><th></th><th colspan="2">Number of IP datapoints</th></tr><tr><th></th><th>Krafla</th><th>Hvedemarken</th></tr></thead><tbody><tr><td>Both rejected</td><td>4109</td><td>1677</td></tr><tr><td>Both accepted</td><td>2334</td><td>2624</td></tr><tr><td>Auto accepted; manual rejected</td><td>864</td><td>434</td></tr><tr><td>Manual accepted; auto rejected</td><td>1008</td><td>1296</td></tr><tr><td>Per cent agree</td><td>77.5 per cent</td><td>71.3 per cent</td></tr><tr><td>Per cent disagree</td><td>22.5 per cent</td><td>28.7 per cent</td></tr><tr><td></td><th colspan="2">Number of IP gates</th></tr><tr><td></td><th>Krafla</th><th>Hvedemarken</th></tr><tr><td>Both removed</td><td>241337</td><td>61945</td></tr><tr><td>Both kept</td><td>35980</td><td>43952</td></tr><tr><td>Discrepancy</td><td>38653</td><td>32816</td></tr><tr><td>Per cent agree</td><td>88 per cent</td><td>76 per cent</td></tr><tr><td>Per cent disagree</td><td>12 per cent</td><td>24 per cent</td></tr></tbody></table>

and a similar peak is observed for Krafla data (Fig. 4g). While Hvedemarken data seem to show a slightly higher value, $\beta_{mean} = 0.4$ (Fig. 4h), the majority of SE models have $\beta$ values falling in the range $\beta = 0.3 – 0.4$.

Similarly to $\alpha$ distribution, the distribution of $\tau_{KWW}$ values show two main peaks. The first peak, observed for all subplots in Fig. 4, corresponds to $\tau_{KWW,mean} = 100$ ms. The second peak, only observed for Nesjavellir data (Figs 4i–j) and mostly for the 64-inch configuration, corresponds to $\tau_{KWW,mean} = 1000$ ms. Most $\tau_{KWW}$ values fall in the range 10–1000 ms for all data sets.

Finally, we stress here that the SE parameters analysed here are not ‘intrinsic’ parameters to the geological media investigated since the SE function used to process the data (eq. 1) corresponds to a theoretical step response and thus does not account for the actual pulse sequence used in these field experiments.

## 4.2 Comparison to manual processing

The manual and automatic processing are compared for Krafla and Hvedemarken data sets. Table 3 shows that there is 70–80
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Spectral inversion of seven DCIP profiles at Krafla, visualized in 3-D for the parameters $\phi_{max}$ (left) and $\tau$ (right) obtained with a Cole-Cole parametrization by the EEMverter inversion software based on the manually (top) and automatically (bottom) processed data sets. Each profile is cut below the depth of investigation. The profile ID is indicated in the corner of the $\phi_{max}$ sections.</figcaption></figure>

per cent agreement for the datapoints and 75-90 per cent agreement for the IP gates. In this section, we further analyse the levels of discrepancy resulting from interpreting the automatically versus manually processed data sets. The details concerning each IP transient that was kept either by the manual or by the automatic processing (negative versus positive gates, gates accepted by the manual versus automatic processing and outcome of the final SE fit) can be inspected in the folders provided as Supplementary Material S1 that contain all the plots, sorted by type of discrepancy following the nomenclature in Table 3.

In general, the automatic processing accepts less datapoints but more gates per datapoint. In Krafla, 1008 out of 8315 datapoints are rejected by the automatic processing while they were accepted by the manual processing; in Hvedemarken, 1296 out of 6031 datapoints are in this case (Table 3). In a minority of these cases, the automatic processing fails to extract the IP signal from the (noisy) discharge curves and this especially happens when there are only a few gates left in the signal accepted by the manual processing, typically less than 10. Yet, in most of these cases, the automatic processing simply rejects IP discharge curves that either show an extreme behaviour among those presented in Fig. 1, such as very flat decays rejected by the criterion $0.2 \le \beta$ (eq. 12) or that correspond to negative decays. On the other hand, the automatic processing also accepts a non-negligible number of datapoints that had been rejected by the manual processing: 864 out of 8315 in Krafla; 434 out of 6031 in Hvedemarken (Table 3). In most cases, these were wrongly rejected by the manual processing, according to a manual inspection of the plots in Supplementary Material S1.

In order to further understand the difference between manual and automatic processing, the seven DCIP profiles in Krafla have been inverted with manual and automatic processing. The inversion results are first compared thanks to a 3-D visualization of the parameters $\phi_{max}$ and $\tau$ obtained at the seven transects with

both processing (Fig. 5). The inverted profiles are generally quite similar with some discrepancy observed. First, it can be observed that the depth of investigation (DOI) is deeper for the automatic processing in the case of ISL10, while it is shallower in the case of ISL1. These two profiles also show the largest discrepancies in terms of parameter spatial distribution: (i) the long relaxation times on the western side of ISL1 with the manual processing are not observed with the automatic processing, (ii) similarly the narrow layer of long relaxation times observed on the northern part of ISL10 with the manual processing is not observed with the automatic processing. High-$\tau$ blocks are observed in both cases in the southern part of ISL10, ISL11 and ISL12, but the contrasts are sharper with automatic processing. Similarly, the high-$\phi_{max}$ layers on these three profiles have a sharper contrast with the automatic processing. In addition, there is a larger block of high-$\phi_{max}$ in the northern part of ISL11 in the case of automatic processing.

While it is impossible to determine which model is better, we also compare for each profile and each processing method the data misfit and the amount of data used in the inversion. A synthesis of this comparison is presented in Table 4, and a more exhaustive comparison is carried out with Figs 6 and 7. Table 4 shows that the data misfit is generally lowered with the automatic processing, except for ISL11 and ISL12 where the automatically processed data set results in a slightly higher misfit, albeit using 15 per cent and 25 per cent more IP gates, respectively, than what remained in the manually processed data sets. It is interesting to note that the inverted transects for ISL11 and ISL12 look fairly similar (Fig. 5). The profile showing the largest discrepancies between the processing methods (ISL10) for the spatial distribution of $\phi_{max}$ and $\tau$, has a lower misfit with the automatic processing (Table 4), and also a larger number of gates that are kept overall (Fig. 7). ISL10 is a profile that contained many negative IP discharge curves, kept during manual
<table><caption>**Table 4.** Comparison of inversion results for the seven DCIP profiles at Krafla, using manual processing and automatic processing. Note that the quadrupoles reported here correspond to a single electrode configuration but in many cases to several IP decays, since different pulse sequences were used for many of the electrode quadrupoles.</caption><thead><tr><th>DCIP profile</th><th>Type of indicator</th><th>Manual Proc.</th><th>Auto Proc.</th></tr></thead><tbody><tr><td rowspan="5">ISL1</td><td># quadrupoles</td><td>824</td><td>823</td></tr><tr><td># IP gates</td><td>4066</td><td>3499</td></tr><tr><td>Data misfit</td><td>4.87</td><td>3.73</td></tr><tr><td>DC misfit</td><td>5.50</td><td>4.27</td></tr><tr><td>IP misfit</td><td>4.70</td><td>3.55</td></tr><tr><td rowspan="5">ISL2</td><td># quadrupoles</td><td>856</td><td>856</td></tr><tr><td># IP gates</td><td>5910</td><td>6349</td></tr><tr><td>Data misfit</td><td>1.80</td><td>1.83</td></tr><tr><td>DC misfit</td><td>2.51</td><td>2.48</td></tr><tr><td>IP misfit</td><td>1.66</td><td>1.701</td></tr><tr><td rowspan="5">ISL3</td><td># quadrupoles</td><td>889</td><td>899</td></tr><tr><td># IP gates</td><td>6001</td><td>6023</td></tr><tr><td>Data misfit</td><td>1.93</td><td>1.72</td></tr><tr><td>DC misfit</td><td>3.23</td><td>3.11</td></tr><tr><td>IP misfit</td><td>1.60</td><td>1.34</td></tr><tr><td rowspan="5">ISL4</td><td># quadrupoles</td><td>928</td><td>928</td></tr><tr><td># IP gates</td><td>7598</td><td>7689</td></tr><tr><td>Data misfit</td><td>1.98</td><td>1.64</td></tr><tr><td>DC misfit</td><td>2.84</td><td>2.82</td></tr><tr><td>IP misfit</td><td>1.83</td><td>1.40</td></tr><tr><td rowspan="5">ISL10</td><td># quadrupoles</td><td>1509</td><td>1509</td></tr><tr><td># IP gates</td><td>12800</td><td>12842</td></tr><tr><td>Data misfit</td><td>3.244</td><td>2.541</td></tr><tr><td>DC misfit</td><td>3.194</td><td>3.153</td></tr><tr><td>IP misfit</td><td>3.251</td><td>2.448</td></tr><tr><td rowspan="5">ISL11</td><td># quadrupoles</td><td>1272</td><td>1272</td></tr><tr><td># IP gates</td><td>10112</td><td>12353</td></tr><tr><td>Data misfit</td><td>2.08</td><td>2.63</td></tr><tr><td>DC misfit</td><td>3.20</td><td>3.33</td></tr><tr><td>IP misfit</td><td>1.86</td><td>2.54</td></tr><tr><td rowspan="5">ISL12</td><td># quadrupoles</td><td>1200</td><td>1200</td></tr><tr><td># IP gates</td><td>9890</td><td>13870</td></tr><tr><td>Data misfit</td><td>2.03</td><td>2.66</td></tr><tr><td>DC misfit</td><td>3.18</td><td>3.77</td></tr><tr><td>IP misfit</td><td>1.81</td><td>2.53</td></tr></tbody></table>

processing and removed during automatic processing, among which a few of them could be fitted with EEMverter inversion software but most of them could not. This data set also contains fast decaying curves that were rejected by the automatic processing—due to the criterion $\beta < 1$ (eq. 12)—and accepted by the manual processing (Fig. 8).

The inversion results in Krafla also allow understanding what happens to the IP decays that were rejected by the manual processing and accepted by the automatic processing. In most cases, they were correctly fitted by the spectral inversion, which illustrates that the manual processing is imperfect. This is not surprising considering that manual processing is usually carried out by inspecting 5 to 10 IP decays together, in order to save time, which may bias the observation of the decays when they are significantly overlapping each other.

Overall, it appears that the automatic processing performs similarly to the manual processing, if not slightly better, while saving significant amount of time. Manual processing of TDIP data is indeed a time-consuming and gruelling task: with an

estimated time of 5 s per datapoint (or IP decay) when carrying out the processing of several datapoints together, it takes 2–3 hr for a typical data set containing 2000 datapoints, provided the user has an appropriate software, which may come as an added cost. On the other hand, automatic processing takes in average 2 s per datapoint (depending on the number of steps necessary) but it is all automatized and can run independently from the user. Manual processing is also subjective and depends on the user, whereas the automatic processing follows clearly formulated and reproducible criteria. Moreover, the automatic processing presented here is a standalone and open code that can be tuned and improved by users, as opposed to automatic processing that may be embedded directly in an inversion software.

### 4.3 Understanding of the integral chargeability behaviour

The behaviour of the integral chargeability for various IP decays is investigated here thanks to the framework offered by SE-modelling of IP decays. Indeed, the SE function allows extending the integration window beyond the field acquisition range, using eq. (7), and in particular calculating an analytical value of the integral assuming an infinite measuring time-window, $[0; +\infty]$ (eq. 15). For clarity purpose in the rest of the article, we refer to the integral chargeability calculated in a time-limited interval as $M_{int,TL}$, and to the analytical integral chargeability $[0; +\infty]$ as $M_{int,+\infty}$.

$$M_{int,+\infty} = \alpha \int_{0}^{+\infty} e^{-(\frac{t}{\tau_{KWW}})^{\beta}} dt = \alpha \frac{\tau_{KWW}}{\beta} \Gamma(\frac{1}{\beta}), \quad (15)$$

where $\Gamma$ is the gamma function (see e.g. E. Artin 2015).

A series of synthetic SE curves, representing a broad range of IP decays, are generated. Their integral chargeability calculated by eqs (7) ($M_{int,TL}$) and (15) ($M_{int,+\infty}$) are compared, considering three different integration intervals for the time-limited calculation (Fig. 9). For the largest integration interval, $[10^{-3} - 10^5]$ s], very limited deviation is observed between $M_{int,TL}$ and $M_{int,+\infty}$, regardless $\beta$ (colours) and $\tau_{KWW}$ (symbol size) values (Fig. 9a). For an intermediate integration interval (Fig. 9b), $[10^{-3} - 10]$ s], which could be considered as a very ambitious target for field measurements, large deviations are found between $M_{int,TL}$ and $M_{int,+\infty}$, unless $\beta$ is close to one. The deviation also worsens with increasing $\tau_{KWW}$ values. For a typical field measurement time-window, $[2 \times 10^{-2} - 2]$ s (Fig. 9c), even larger deviations occur; deviations become larger with decreasing $\beta$. Note that, for $\beta$ close to one, deviations are still observed for $\tau_{KWW}$ values outside the integration range.

It is evident from the synthetic analysis that the integration window, which is normally determined based on the measuring window for the voltage signal, plays an essential role on the time-limited integral value (Fig. 9). The discrepancy between $M_{int,TL}$ and $M_{int,+\infty}$ becomes larger as the integration window decreases and as $\beta$ decreases. Unless the integration window is unrealistically large, the values are only equivalent in cases where $\beta$ is close to one, and $\tau_{KWW}$ is within the integration time range, that is, for a limited range of $\tau_{KWW}$. Since field measurements impose a sharp limit to the acquisition time-range, the question whether field data can be artificially extended beyond the time-window used in the field—thanks to the SE fitted function, is legitimate. Artificially extending the intervals thanks to the SE
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Pseudo-sections of the seven DCIP profiles at Krafla, showing the IP data misfit after spectral inversion for the manual (left) and automatic (right) processing.</figcaption></figure>

function would allow two things in particular: (i) calculating a time-limited integral chargeability over consistent integration windows regardless the actual data range (which depends on the current waveform, the duty cycle, as well as the gating strategy, the voltage sampling rate and earliest and latest gates left after the processing) and (ii) calculating an RTD over a larger range of relaxation times. Therefore, we explore here to what extent this approach is susceptible to bias the analysis, considering that (i) the link between an SE-model and a field IP decay is non-unique, both due to data noise and to the limited time-range of the field data and (ii) the simple SE-fitting approach proposed here does not take into account the current waveform.

We carry out this exploration using the Nesjavellir field data because (i) a larger number of datapoints are available and (ii) the logging configuration allows evaluating the consistency over depth, while it is more challenging with cross-borehole or surface configurations. The SE functions fitted on Nesjavellir data are used for calculating both $M_{int,+\infty}$ (eq. 15) and $M_{int,TL}$ (eq. 7),

similarly to Fig. 9. For the latter, this means that polarizability values $\eta(t)$ (eq. 8) are recalculated based on the fitted SE function at each gate, instead of using the field polarizability data. This allows both to smoothen out noise and also to test the behaviour of $M_{int,TL}$ when gradually extending the integration window beyond the field acquisition window.

As shown in Fig. 10, a much better consistency over depth is observed with $M_{int,TL}$ calculated on the field data-range [1.1 ms; 2 s], than with $M_{int,+\infty}$. While the consistency between $M_{int,TL}$ and $M_{int,+\infty}$ improves when gradually extending the integration window for $M_{int,TL}$ (in agreement with results presented in Fig. 9), the scattering of $M_{int,TL}$ as a function of depth worsens, leaning towards the scattering observed for $M_{int,+\infty}$ (Fig. 10). The latter observation could not have been predicted from the synthetic analysis itself.

The increased depth-scattering suggests that the SE-fitted parameters are underdetermined for IP data acquired on a limited time-range, considering that these data also carry noise. In other
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Pseudo-sections of the seven DCIP profiles at Krafla, showing the number of used gates per quadrupole in the spectral inversion for the manual (left) and automatic (right) processing. The grey dots indicate quadrupoles included in the DC inversion but not the spectral IP inversion (not enough gates left after processing).</figcaption></figure>

words, the SE models are non-unique and the artificial extension of the integration window beyond the original data time-range may introduce a significant bias. Therefore, using the SE functions beyond their outlier-removal processing role may not be appropriate. From this combined synthetic and field data analysis, we conclude that (i) the SE functions should not be used outside the field acquisition time-range and (ii) the theoretical $M_{int,+\infty}$ value should not be used for interpreting field data.

## 4.4 Using DD parameters for interpreting spectral IP changes in Nesjavellir data set

Debye decomposition is carried out on Nesjavellir data set, using directly the SE fitted functions, which allows minimizing the influence of data noise on the DD parameters; this effect is

described in details in L. Lévy *et al.* (2024, Appendix B, Figs. B1 and B2). The data time-range chosen for DD is [3.3; 1962 ms] and IP decays that are not reasonably fitted within this interval are not processed further.

The relaxation time interval for DD is fixed accordingly to [3.3; 1962 ms], which corresponds in practice to a set of $N = 300$ logarithmically-spaced values of discrete relaxation times $\tau_{k,k=1..N}$ in the range 3.3 – 1962 ms. IP decays fitted by a SE model with $\tau_{KWW}$ outside this interval are discarded. When considering the whole Nesjavellir data set, that is including wells NN3, NN4, NN7, as well as both 16- and 64-inch spacing, 90 per cent of IP decays meet the two conditions (well fitted by the SE function in the chosen interval, and $\tau_{KWW}$ within this interval) out of the 24187 IP decays that enter step #6 of the processing algorithm (Table 2).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Negative transients kept during the manual processing of transect ISL10: IP discharge curves showing negative data and forward simulations together with the corresponding quadrupole positions on the IP data misfit pseudosection. The red colour indicates that the data are negative (similarly to Fig. 3)—their absolute value is shown here. Each subplots shows data for two different quadrupoles sharing the same pseudoposition and there is sometime more than one datapoint per quadrupole (if different pule sequences are used for the same quadrupole).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of the integral chargeability value, when calculated on a time-limited interval, $M_{int,TL}$ and with the analytical formula, $M_{int,+\infty}$, for a broad range of IP decays and of integrating windows. The value of $\alpha = 500\text{mV V}^{-1}$ is fixed. The colour scale corresponds to $\beta$ varying from 0.2 to 1 and the symbol size to $\tau_{KWW}$ varying from 50 to 5000 ms. The three sub-figures correspond to three integration windows, from left to right: $[10^{-3} - 10^5 \text{s}]$, $[10^{-3} - 10 \text{s}]$ and $[2 \times 10^{-2} - 2 \text{s}]$.</figcaption></figure>

In order to exploit the spectral information embedded in the DD, the relaxation time interval is divided into two sub-intervals, 3.3–100 ms and 100–1962 ms, and the $m_k$ coefficients (eq. 13) in these two intervals are added to obtain $m_{low}$ and $m_{high}$, respectively. Fig. 11(b) shows that $m_{low}$ does not change over time in NN4, while $m_{high}$ shows significant time-variations at +40 d. The spatial variations are consistent over time for both $m_{low}$ and $m_{high}$. Similar observation can be made for NN3 (Fig. 11a), especially below 450 m depth. It is interesting to note that, in NN3, the section below 450 m has relatively small $m_{low}$ values at all dates compared to the section above, while the $m_{high}$ values become larger below 450 m than in the rest of the well at +40 d.

These qualitative observations may indicate that, in NN4 as well as below 450 m in NN3, the volume of particles that polarize at longer time (i.e. larger particles) increases over time, while the volume of particles polarizing at shorter times (smaller particles) remains more or less constant. The integral chargeability alone would not allow such interpretations. However, further quantitative interpretation of the spectral IP changes in terms of petrophysical and mineralogical parameters would require convolving the step response with the actual current waveform.

## 4.5 Comparing DD and SE parameters

First, it is worthwhile to note that, for $\beta = 1$ the SE function presented in eq. (1) describes an exponential decay corresponding to one Debye element, which can also be represented by a RC circuit (see eq. 13). Then, Fig. 12 shows that the ratios $\frac{m_{tot}}{\alpha}$ and $\frac{\tau_{RTD,mean}}{\tau_{KWW}}$ are close to one for $\beta > 0.5$ and $\beta > 0.4$, respectively (Figs 12a–c and d–f). For $\beta < 0.4 - 0.5$, we see that $\frac{m_{tot}}{\alpha} < 1$ while $\frac{\tau_{RTD,mean}}{\tau_{KWW}} > 1$. In other words, if the SE model parameters $\alpha$ and $\tau_{KWW}$ were directly taken in place of the DD mean parameters in case of $\beta < 0.5$, they would overestimate $m_{tot}$ and underestimate of $\tau_{RTD,mean}$, respectively. On the other hand, a fairly reliable prediction of $m_{tot}$ can be obtained using both $\alpha$ and $\beta$ through an empirical function of the form $m_{tot} = \alpha \text{erf}(a \cdot \beta - b)$, fitted to the trends observed on Figs 12(a)–(c).

## 5 DISCUSSION: LIMITATIONS OF THE AUTOMATIC PROCESSING PROCEDURE

A choice is made in the processing to ignore gates with negative chargeability data. Therefore, the approach is not suitable for data sets containing significant number of datapoints with
NN3 - 64”

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Comparison of $M_{\text{int,TL}}$ (calculated on the field data-range, that is, [1.1 ms; 2 s]) and $M_{\text{int},+\infty}$ depth-consistency for Nesjavellir data in NN3 wells with 64-inch electrode spacing at the four monitoring dates (four colours).</figcaption></figure>

negative chargeability decays or other types of heterodox transients where negative gates play an important role (G. Fiandaca et al. 2022), provided that these negative data can actually be fitted by an inversion software. Nevertheless, we show here that even state-of-the-art spectral IP inversion software struggle to fit these heterodox transients (Figs 6 and 8).

More generally, the KWW model (referred to here as the SE function) covers a large range of phenomenological behaviours—larger than the Cole-Cole model—but does not cover all possible shapes of TDIP signals. In addition, decays with low $\beta$ are uncertain to fit on realistic field time-ranges, which has led to ignoring datapoints fitted by an SE model with $\beta < 0.2$. The automatic processing does reject heterodox transients that could in theory be fitted by spectral inversion (G.

Fiandaca et al. 2022). However, we show that, in the Krafla data set, negative decays are not well fitted in practice (Fig. 8) and rather contribute to increasing the data misfit without improving the resolution. On the other hand, the automatic processing accepts a fair amount of IP decays that were wrongly rejected by the manual processing. The automatic processing also tends to keep more gates per IP datapoint that are correctly fitted during the inversion (see Table 4 and Fig. 7 for the Krafla data sets). The detail of all the cases of discrepancy between the manual and automatic processing can be inspected in the folders provided as Supplementary Material S1 and containing all the corresponding plots.

An empirical detection method for the characteristic EM 'kink' point is proposed here in order to handle the EM noise, which in some cases also corresponds to the effect of low-pass digital filters. This step, which is solely based on the user experience during manual processing, significantly improves the performance of the SE-based processing procedure. Nevertheless, further investigation of this characteristic point should be warranted to better understand EM coupling effects in TDIP and in particular to evaluate the amplitude of remaining EM coupling effects after this kink point.

All this being considered, the overwhelming advantage of the automatic processing is that it is fully automatic (saving time to the user), with transparent and objective criteria (as opposed to the subjectivity of a manual user) and it is a standalone procedure that does not require an inversion software attached to it. In the future, case-specific improvements may allow better handling of the few failures of the automatic processing. This may include, for example, for TDIP wireline logging data, the averaging of IP decays over a certain depth range (see e.g. L. Lévy et al. 2024), as well as starting the SE fit using the fitting result of a neighbouring depth as the initial set of parameters.

## 6 CONCLUSION

A new TDIP data processing method is proposed and the corresponding MATLAB scripts are shared as open-source codes. The full procedure, based on the KWW function (R. Kohlrausch 1854; G. Williams & D.C. Watts 1970), as well as an empirical handling of the early time electromagnetic coupling effects, is tested and validated on three data sets that cover a large range of contexts, electrode configurations and acquisition settings. These three data sets include logging, cross-borehole and surface TDIP configurations, with electrode spacing ranging from a few centimetres to several tens of metres. The procedure allows smoothing out the TDIP decays as well as excluding outlier decays, while keeping a large number of them. Compared to simply calculating the integral chargeability, this approach is especially useful for further extracting spectral information from the TDIP data, while keeping it simple for unexperienced users that only have a gated IP signal provided by the instrument at their disposal. Indeed, while DD is often presented as the most suitable approach for converting TDIP data into meaningful spectral parameters, this approach is very sensitive to data noise and can easily act as a noise amplifier. The procedure proposed here allows overcoming this problem in an automatic and transparent way by first processing the data, before applying a Debye decomposition. Moreover, a clear and consistent correlation is observed between the SE parameter $\alpha$ and the DD parameter $m_{\text{tot}}$ across
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Depth and time variations of $m_{low}$ and $m_{high}$ for NN3 and NN4.</figcaption></figure>

the three data sets. Similarly, a correlation is observed between $\tau_{KWW}$ and $\tau_{mean,RTD}$.

The KWW function also offers a theoretical framework for exploring (i) the integral chargeability behaviour when varying the integration window and the decay shape and (ii) the information embedded in the SE fitting parameters. It is clear, from the Nesjavellir logging TDIP data, that the integral chargeability calculated on the field data time-range is a relevant qualitative parameter that allows efficiently smoothing out the noise, as illustrated by the consistent variations in depth and time. The progressive extension of the data time-range, using the SE-fitted function, beyond the original acquisition range leads to increasing depth-scattering. Similar depth-scattering is observed with the infinite integral chargeability value, $M_{int,+\infty}$, calculated with the analytical integral formula of the KWW function. These field-based observations highlight the fact that the SE-fitted function may be non-unique, due to the limited time-range available for fitting and the presence of noise, and suggest that the SE functions should not be used beyond their outlier-removal processing role. Nevertheless, the SE functions allow investigating the influence of data time-range and relaxation time-range on DD parameters through synthetic analysis of SE-based TDIP decays. We could thus propose recommendations on how to tune the different

settings when carrying out a DD on IP decays, for a given data set. Finally, we show that DD-extracted spectral information can be used beyond the $m_{low}$ and $m_{high}$, as they allow further interpretation in terms of mineralization changes over time. In addition, the SE and DD parameters obtained here without convolution with the current waveform are not intrinsic to the investigated medium cannot be used for quantitative conversion to petrophysical parameters. Nevertheless, we propose to use the spectral IP information extracted from the DD in a qualitative manner, through the parameters $m_{low}$ and $m_{high}$, to evaluate whether IP changes happening over time in case of time-lapse wireline logging TDIP data (Nesjavellir data set) correspond to IP changes at early or late relaxation times.

## DATA AVAILABILITY

The open source code used for data processing, as well as the input data and output parameters can be found on the Open Access Repository 'EaSy Data', through the following link: https://doi.org/10.57932/f28974f0-5c2c-4731-9925-0c5145b9998d.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Relation between the ratio $\frac{m_{tot}}{\alpha}$ and $\beta$ (a–c) and between the ratio $\frac{\tau_{RTD,mean}}{\tau_{KWW}}$ and $\beta$ (d–f) for the three data sets.</figcaption></figure>

## SUPPORTING INFORMATION

Supplementary data are available at *GJIRAS* online.

**Supplementary Material S1:** zip folder containing a series of .png plots sorted into six sub-folders, three for Hvedemarken data set and three for Krafla data set, in each case corresponding to (i) IP decays accepted by manual and automatic processing, (ii) IP decays accepted by automatic processing and rejected by manual processing, (iii) IP decays accepted by manual processing and rejected by automatic processing.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## ACKNOWLEDGMENTS

LL thanks Gianluca Fiandaca and Nicole Sullivan for their help using the latest version of EEMverter inversion software, as well as EEMstudio visualization tool. The authors warmly thank Dr Seogi Kang and another anonymous reviewer for their thorough feedback and thus significant contribution to the final manuscript.

## REFERENCES

Alvarez, F., Alegra, A. & Colmenero, J., 1991. Relationship between the time-domain Kohlrausch-Williams-Watts and frequency-domain Havriliak-Negami relaxation functions, *Phys. Rev. B*, **44**, 7306–7312.

Artin, E., 2015. *The Gamma Function*, Courier Dover Publications.

Barfod, A.S., Lévy, L. & Larsen, J.J., 2021. Automatic processing of time domain induced polarization data using supervised artificial neural networks, *Geophys. J. Int.*, **224**(1), 312–325.

Bertin, J. & Loeb, J., 1976. *Experimental and Theoretical Aspects of Induced Polarization*, Vol. 1, Schweizerbart Science Publishers.

Binley, A. & Slater, L., 2020. *Resistivity and Induced Polarization: Theory and Applications to the Near-Surface Earth*, Cambridge Univ. Press.

Cole, K.S. & Cole, R.H., 1941. Dispersion and absorption in dielectrics i. alternating current characteristics, *J. Chem. Phys.*, **9**(4), 341–351.

Cox, L.H., Zhdanov, M.S. & Alfouzan, F., 2019. Inductive coupling in IP measurements and applications to 3D distributed array field data inversion, in *SEG Technical Program Expanded Abstracts 2019*, pp. 2173–2177, Society of Exploration Geophysicists. doi: 10.1190/segam2019-3213081.1.

Dahlin, T. & Leroux, V., 2012. Improvement in time-domain induced polarization data quality with multi-electrode systems by separating current and potential cables, *Near Surf. Geophys.*, **10**(6), 545–565.

Davidson, D.W. & Cole, R.H., 1951. Dielectric relaxation in glycerol, propylene glycol, and n-propanol, *J. Chem. Phys.*, **19**(12), 1484–1490.
Dey, A. & Morrison, H.F., 1973. Electromagnetic coupling in frequency and time-domain induced-polarization surveys over a multilayered earth, *Geophysics*, **38**(2), 380-405.

Fiandaca, G., Auken, E., Christiansen, A.V. & Gazoty, A., 2012. Time-domain-induced polarization: full-decay forward modeling and 1D laterally constrained inversion of cole-cole parameters, *Geophysics*, **77**(3), E213-E225.

Fiandaca, G., Ramm, J., Binley, A., Gazoty, A., Christiansen, A.V. & Auken, E., 2013. Resolving spectral information from time domain induced polarization data through 2-D inversion, *Geophys. J. Int.*, **192**, 631-646.

Fiandaca, G., Olsson, P.I., Maurya, P.K., Kühl, A., Bording, T., Dahlin, T. & Auken, E., 2022. Heterodox transients in time-domain-induced polarization, *Geophysics*, **87**(1), E35-E47.

Fiandaca, G., Chen, J. & Zhang, B., 2024. Eemverter, a new modelling tool for electric and electromagnetic data with focus on induced polarization, in *NSG 2024 30th European Meeting of Environmental and Engineering Geophysics*, vol. **2024**, pp. 1-5, European Association of Geoscientists and Engineers, Helsinki, Finland.

Flores Orozco, A., Williams, K.H., Long, P.E., Hubbard, S.S. & Kemna, A., 2011. Using complex resistivity imaging to infer biogeochemical processes associated with bioremediation of an uranium-contaminated aquifer, *J. geophys. Res.: Biogeosciences*, **116**(G3). doi: 10.1029/2010JG001591.

Flores Orozco, A., Gallistl, J., Bücker, M. & Williams, K.H., 2018. Decay curve analysis for data error quantification in time-domain induced polarization imaging, *Geophysics*, **83**(2), E75-E86.

Fullagar, P.K., Zhou, B. & Bourne, B., 2000. EM-coupling removal from time-domain IP data, *Explor. Geophys.*, **31**(1-2), 134-139.

Gazoty, A., Fiandaca, G., Pedersen, J., Auken, E. & Christiansen, A.V., 2013. Data repeatability and acquisition techniques for time-domain spectral induced polarization, *Near Surf. Geophys.*, **11**(4), 391-406.

Havriliak, S. & Negami, S., 1966. A complex plane analysis of α-dispersions in some polymer systems, *J. Polymer Sci.: Polymer Symp.*, **14**, 99-117.

He, J.S., Xiong, B., BAO, L.Z. & YU, Y.C., 2008. Time characteristics of EM coupling in induced polarization, *Chin. J. Geophys.*, **51**(3), 645-653.

Hilfer, R., 2002. Part 1-structured and complex fluids-h-function representations for stretched exponential relaxation and non-debye susceptibilities in glassy systems, *Phys. Rev.-Section E-Stat. Nonlin. Soft Matter Phys.*, **65**(6), 61 510-61510.

Holm, S., 2020. Time domain characterization of the cole-cole dielectric model, *J. Electric. Bioimpedance*, **11**(1), 101.

Izumoto, S. et al., 2022. Pore-scale mechanisms for spectral induced polarization of calcite precipitation inferred from geo-electrical millifluidics, *Environ. Sci. Technol.*, **56**(8), 4998-5008.

Kang, S. & Oldenburg, D.W., 2019. Inversions of time-domain spectral induced polarization data using stretched exponential, *Geophys. J. Int.*, **219**(3), 1851-1865.

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

Kessouri, P. et al., 2019. Induced polarization applied to biogeophysics: recent advances and future prospects, *Near Surf. Geophys.*, **17**(6), 595-621.

Kohlrausch, R., 1854. Theorie des elektrischen Rückstandes in der Leidener Flasche, *Ann. Phys.*, **167**(2), 179-214.

Lévy, L. et al., 2019. Electrical resistivity tomography and time-domain induced polarization field investigations of geothermal areas at Krafla, Iceland: comparison to borehole and laboratory frequency-domain electrical observations, *Geophys. J. Int.*, **218**(3), 1469-1489.

Lévy, L. et al., 2022. Quantifying reagent spreading by cross-borehole electrical tomography to assess performance of groundwater remediation, *Water Resour. Res.*, **58**(9), e2022WR032218..

Lévy, L., Ciraula, D.A., Legros, B., Martin, T. & Weller, A., 2024. Understanding the fate of H₂S injected in basalts by means of time-domain induced polarization geophysical logging, *J. geophys. Res.: Solid Earth*, **129**(6), e2023JB028316..

Li, Y. & Oldenburg, D.W., 1999. 3-D inversion of DC resistivity data using an l-curve criterion, in *SEG International Exposition and Annual Meeting*, SEG-1999, SEG. doi: 10.1190/1.1820992.

Mao, D., Revil, A. & Hinton, J., 2016. Induced polarization response of porous media with metallic particles part 4: detection of metallic and nonmetallic targets in time-domain-induced polarization tomography, *Geophysics*, **81**(4), D359-D375.

Martin, T., Günther, T., Orozco, A.F. & Dahlin, T., 2020. Evaluation of spectral induced polarization field measurements in time and frequency domain, *J. Appl. Geophys.*, **180**, 104141.

Martin, T., Titov, K., Tarasov, A. & Weller, A., 2021. Spectral induced polarization: frequency domain versus time domain laboratory data, *Geophys. J. Int.*, **225**(3), 1982-2000.

Maurya, P.K., Balbarini, N., Møller, I., Ronde, V., Christiansen, A.V., Bjerg, P.L., Auken, E. & Fiandaca, G., 2018. Subsurface imaging of water electrical conductivity, hydraulic permeability and lithology at contaminated sites by induced polarization, *Geophys. J. Int.*, **213**(2), 770-785.

Mellage, A., Smeaton, C.M., Furman, A., Atekwana, E.A., Rezanezhad, F. & Van Cappellen, P., 2018. Linking spectral induced polarization (SIP) and subsurface microbial processes: results from sand column incubation experiments, *Environ. Sci. Technol.*, **52**(4), 2081-2090.

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

Olsson, P.I., Dahlin, T., Fiandaca, G. & Auken, E., 2015. Measuring time-domain spectral induced polarization in the on-time: decreasing acquisition time and increasing signal-to-noise ratio, *J. Appl. Geophys.*, **123**, 316-321.

Olsson, P.I., Fiandaca, G., Larsen, J.J., Dahlin, T. & Auken, E., 2016. Doubling the spectrum of time-domain induced polarization by harmonic de-noising, drift correction, spike removal, tapered gating and data uncertainty estimation, *Geophys. J. Int.*, **207**(2), 774-784.

Pelton, W., Sill, W. & Smith, B., 1983. Interpretation of complex resistivity and dielectric data. Part I, *Geophys. Trans.*, **29**(4), 297-330.

Pelton, W.H., Ward, S.H., Hallof, P.G., Sill, W.R. & Nelson, P.H., 1978. Mineral discrimination and removal of inductive coupling with multifrequency IP, *Geophysics*, **43**(3), 588-609.

Radic, T., 2004. Elimination of cable effects while multi-channel sip measurements, in *Near Surface 2004-10th EAGE European Meeting of Environmental and Engineering Geophysics*, Utrecht, Netherlands, pp. cp-10, European Association of Geoscientists and Engineers.

Rembert, F., Stolz, A., Soulaine, C. & Roman, S., 2023. A microfluidic chip for geoelectrical monitoring of critical zone processes, *Lab on a Chip*, **23**(15), 3433-3442.

Schlumberger, C., 1920. Etude sur la prospection électrique du sous-sol, Gauthier-Villars (Paris), pp 1-115.

Sumner, J.S., 1976. *Principles of Induced Polarization for Geophysical Exploration*, Elsevier Scientific Publishing Company Amsterdam.

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

Tarasov, A. & Titov, K., 2013. On the use of the cole-cole equations in spectral induced polarization, *Geophys. J. Int.*, **195**(1), 352-356.

Telford, W.M., Geldart, L.P. & Sheriff, R.E., 1990. *Applied Geophysics*, 2nd edn, Cambridge Univ. Press.

Wang, H., Huisman, J.A., Zimmermann, E. & Vereecken, H., 2021. Experimental design to reduce inductive coupling in spectral electrical impedance tomography (SEIT) measurements, *Geophys. J. Int.*, **225**(1), 222-235.

Wang, H., Huisman, J.A., Zimmermann, E. & Vereecken, H., 2024. Tackling capacitive coupling in broad-band spectral electrical impedance tomography (SEIT) measurements by selecting electrode configurations, *Geophys. J. Int.*, **238**(1), 187-198.

Weller, A. & Slater, L., 2022. Ambiguity in induced polarization time constants and the advantage of the Pelton model, *Geophysics*, **87**(6), E393-E399.
Williams, G. & Watts, D.C., 1970. Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function, *Trans. Farad. Soc.*, **66**, 80–85.

Williams, K.H. et al., 2009. Geophysical monitoring of coupled microbial and geochemical processes during stimulated subsurface bioremediation, *Environ. Sci. Technol.*, **43**(17), 6717–6723.

Zimmermann, E., Kemna, A., Berwix, J., Glaas, W., Mnch, H.M. & Huisman, J.A., 2008. A high-accuracy impedance spectrometer for measuring sediments with low polarizability, *Meas. Sci. Technol.*, **19**(10), 105603.