Relationship between Cole–Cole model parameters and spectral decomposition parameters derived from SIP data

SUMMARY Spectral induced polarization (SIP) data are commonly analysed using phenomenological models. Among these models the Cole–Cole (CC) model is the most popular choice to describe the strength and frequency dependence of distinct polarization peaks in the data. More ﬂexibility regarding the shape of the spectrum is provided by decomposition schemes. Here the spectral response is decomposed into individual responses of a chosen elementary relaxation model, mathematically acting as kernel in the involved integral, based on a broad range of relaxation times. A frequently used kernel function is the Debye model, but also the CC model with some other a priorly speciﬁed frequency dispersion (e.g. Warburg model) has been proposed as kernel in the decomposition. The different decomposition approaches in use, also including conductivity and resistivity formulations, pose the question to which degree the integral spectral parameters typically derived from the obtained relaxation time distribution are biased by the approach itself. Based on synthetic SIP data sampled from an ideal CC response, we here investigate how the two most important integral output parameters deviate from the corresponding CC input parameters. We ﬁnd that the total chargeability may be underestimated by up to 80 per cent and the mean relaxation time may be off by up to three orders of magnitude relative to the original values, depending on the frequency dispersion of the analysed spectrum and the proximity of its peak to the frequency range limits considered in the decomposition. We conclude that a quantitative comparison of SIP parameters across different studies, or the adoption of parameter relationships from other studies, for example when transferring laboratory results to the ﬁeld, is only possible on the basis of a consistent spectral analysis procedure. This is particularly important when comparing effective CC parameters with spectral parameters derived from decomposition results.


I N T RO D U C T I O N
Over the last decade the spectral induced polarization (SIP) method has seen a rapid increase of its use in hydrogeological and environmental studies (e.g.Kemna et al. 2012).SIP data are given as frequency dependent electrical impedance or admittance measurements that can be converted to resistivities or conductivities by accounting for the measurement geometry.
Two basic types of phenomenological models are in use for the quantitative description of SIP data.The Cole-Cole (CC) model and its variants (e.g.Cole & Cole 1941;Pelton et al. 1978;Dias 2000) are typically used if the spectrum exhibits a distinct peak.Here three parameters, the CC relaxation time, the CC chargeability and the CC frequency exponent, account for the position of the peak along the frequency axis, the magnitude and the frequency dispersion of the response, respectively.
Alternatively, SIP data can be decomposed into individual responses of a chosen elementary relaxation model, mathematically acting as kernel in the involved integral, based on a broad range of relaxation times.A frequently used kernel function is the Debye model (e.g.Uhlmann & Hakim 1971;Lesmes & Morgan 2001;Tarasov et al. 2003;Nordsiek & Weller 2008), the special case of the CC model with strongest possible frequency dispersion.However, also the CC model with some other a priorly specified frequency dispersion (e.g. the Warburg model with intermediate frequency dispersion) has been proposed as kernel in the decomposition (Tarasov et al. 2003;Florsch et al. 2014;Revil et al. 2014).The choice of the latter follows pore-scale modelling results, which indicate that the elementary polarization response of a single pore or grain exhibits a broader frequency dispersion than the Debye model response (e.g.Wong 1979;Titov et al. 2002;Bücker & Hördt 2013).Therefore a CC decomposition (CCD) (with fixed frequency dispersion) seems to be more adequate from a petrophysical point of view.The decomposition approach yields a relaxation time distribution (RTD) from which integral spectral parameters similar to those of the CC model can be computed, especially a mean relaxation time and a total chargeability.
It is important to note that the decomposition of SIP spectra based on CC-type kernel functions is highly ambiguous and the resultant RTD depends on the frequency dispersion of the kernel.For instance the equivalent representation of a CC model response as a superposition of Debye model responses based on an appropriate RTD was already recognized by Cole & Cole (1941).However, even if the frequency dispersion of the kernel is fixed in the decomposition, which is commonly done due to this inherent equivalence, the resulting inverse problem typically still is ill-posed, requiring regularization (Florsch et al. 2012(Florsch et al. , 2014;;Weigand & Kemna 2016).This leads to considerable uncertainties in the decomposition results on top of the contribution on the account of data errors.Several studies have assessed uncertainties in the estimation of spectral parameters for the CC model (Ghorbani et al. 2007;Chen et al. 2008) and the Debye decomposition (DD; Florsch et al. 2012;Keery et al. 2012) in detail.
With the growing number and sophistication of SIP studies relating phenomenological model parameters to petrophysical properties, for example to permeability (e.g.Binley et al. 2005;Revil & Florsch 2010;Zisser et al. 2010;Weller et al. 2015), consistency of data analysis procedures becomes important to ensure comparability of results.However, the variety of different models in use today is quite large, and formulations based on either resistivity or conductivity further complicate the interpretation of results across studies.For example, Revil et al. (2015) used four different approaches to determine characteristic relaxation times from SIP data.As pointed out by Weller et al. (2015), the total chargeability derived from the DD is biased by the shape of the spectrum and the analysed frequency range, while the CC chargeability accounts for all frequencies.Gurin et al. (2015) also mention that the total chargeability from the DD is not comparable to the CC chargeability, and is dependent on the position of the spectral peak and the frequency range covered by the data.
Despite these insights, to our knowledge it has not yet been systematically studied how the total chargeability and the mean relaxation time derived from the CC decomposition (including the Debye and Warburg cases) of CC-type data deviate from the original CC chargeability and CC relaxation time.We here investigate these relationships based on synthetic SIP data sampled from an ideal CC response for a broad range of CC parameter values.The results of our study yield quantitative information on the possible over-or underestimation of chargeability and relaxation time when different models are used to analyse SIP data.

The Cole-Cole model
Based on the original formulation of Cole & Cole (1941), the complex conductivity ( σ ) formulation of the CC model is given as (e.g. with ω denoting the angular frequency, σ 0 and σ ∞ the conductivity in the low-and high-frequency limit, respectively, m = (σ ∞ − σ 0 )/σ ∞ the chargeability, j the imaginary unit, τ (σ ) the CC relaxation time, and c the CC frequency exponent.An alternative formulation in terms of complex resistivity ( ρ) was proposed by Pelton et al. (1978) as with ρ 0 = 1/σ 0 being the DC resistivity, ρ ∞ = 1/σ ∞ , and τ (ρ) the CC relaxation time in the resistivity formulation.Importantly, from ρ(ω) = 1/ σ (ω) it follows that τ (σ ) and τ (ρ) are not identical but related by (Florsch et al. 2012) It is easy to show that the peak frequencies of the imaginary components of σ (ω) and − ρ(ω) are related to the relaxation times τ (σ ) and τ (ρ) according to and f . (5) In combination with eq. ( 4) this implies that σ (ω) and −ρ (ω) peak at different frequencies (with denoting the imaginary component).

The Cole-Cole decomposition
Based on a broad range of relaxation times spanned by N discretized values τ k , the discrete form of the CCD in conductivity and resistivity formulation, respectively, can be written as where the kernel (assumed elementary relaxation model) is adopted from the CC model in eqs ( 1) and ( 2), respectively.In eqs ( 6) and ( 7), m k denote the chargeability weights at the sampled relaxation times τ k in conductivity and resistivity formulation, respectively.Usually the relaxation time range is chosen to cover at least the measurement frequency (f) range of the data according to the inverse relationship τ = 1/(2π f) (e.g.Weigand & Kemna 2016).The frequency dispersion of the kernel functions in the decompositions is controlled by the chosen fixed value for c, with the DD resulting for c = 1 and a Warburg decomposition (Revil et al. 2014) resulting for c = 0.5.
The RTD is given by the distribution m k (τ k ), respectively, from which the following integral parameters can be derived: (i) The total chargeability m tot = N k=1 m k (e.g.Nordsiek & Weller 2008) is the analogue to the chargeability m in eqs (1) and ( 2) when a single CC model is used to describe SIP data.However, in the typical case that the fixed value of c in the decomposition is considerably larger than the c value of a single CC model describing the same SIP data, that is, if the elementary relaxation model has a much stronger frequency dispersion than the given data, m tot primarily accounts for polarization within the considered frequency range of the decomposition (spanned by the chosen τ k values), whereas m also contains significant contributions from outside this range.
(ii) The median relaxation time τ 50 is the relaxation time at which 50 per cent of the total chargeability is reached (Nordsiek & Weller 2008;Zisser et al. 2010).
(iii) The mean logarithmic relaxation time τ mean is the chargeability-weighted logarithmic mean value of the RTD (Nordsiek & Weller 2008): ).
(iv) The arithmetic mean of the RTD is given by τ a = N k=1 m k τ k N k=1 m k (Tong et al. 2004).
Although originally proposed for the DD, the above parameters can be computed analogously from the RTD of the CCD.

N U M E R I C A L E X P E R I M E N T S
Synthetic SIP data were generated assuming a CC model response in resistivity formulation (eq.2) by systematically varying τ (ρ) in the range from 1.59 × 10 −5 s to 159.15 s and c from 0.05 to 1. Assumed measurement frequencies were selected using eq.( 5) such that the polarization peaks lie within the frequency range covered by the data, with 30 frequencies in total.CCDs of the resulting spectra were performed for the c values 0.3, 0.5 and 1.0, according to eq. ( 7) (resistivity formulation), using the open-source code from Weigand & Kemna (2016), with 20 relaxation times per decade, and two decades of relaxation times outside the low-and high-frequency limits of the data (222 relaxation times in total, 140 representing the data frequency range).
The input CC chargeability m was compared to the m tot parameters computed from the CCD results.For the relaxation times, the three parameters τ mean , τ 50 , and τ a based on the decomposition with c = 1 (DD) were related to the input CC parameter τ (ρ) .In addition, τ mean resulting from the CCD for the c values 0.3 and 0.5, respectively, was compared to τ (ρ) .Finally, the analogue of eq. ( 4) for the mean relaxation times τ (σ )  mean and τ (ρ) mean obtained from the DD (c = 1) in conductivity (eq.6) and resistivity (eq.7) formulation, respectively, was evaluated for different input CC chargeability values in terms of τ (σ )  mean /τ (ρ) mean .

Chargeability recovery
Original CC chargeabilities m are increasingly underestimated by m (ρ) tot by up to 80 per cent for decreasing c values, that is, decreasing frequency dispersion of the data, and for τ (ρ) approaching either the low or high relaxation time limit in the decomposition (Fig. 1).This behaviour is found for different c values in the range c < c.

Relaxation time recovery
The relaxation times τ (ρ)  mean and τ (ρ) 50 recovered from the CCD (eq.7) over-or underestimate the original CC relaxation time τ (ρ) by up to three orders of magnitude for τ (ρ) approaching the high or low relaxation time limit, respectively, and for decreasing c values (Fig. 2).The absolute deviation patterns are similar to those obtained for m (ρ) tot (Fig. 1).The arithmetic mean relaxation time (τ (ρ)  a ) shows only a narrow band of CC parameters where the agreement with the original τ (ρ) value is good; outside this range, deviations up to three orders of magnitude exist.
The deviation of the CCD-derived (eq.7) relaxation time τ (ρ)  mean for different values of c exhibits a similar pattern as the corresponding m (ρ) tot results (Fig. 1).The absolute deviation from the original τ (ρ) value increases with decreasing c (< c) and with τ (ρ) approaching one of the relaxation time limits (Fig. 3).
The ratio between the DD-derived relaxation times τ (σ ) mean and τ (ρ) mean (using eqs 6 and 7, respectively) diverts from the behaviour of the CC counterparts τ (σ ) and τ (ρ) (described by eq. 4) and shows much less dependence on m (Fig. 4).However, the ratio deviates from 1 (the value which indicates perfect agreement between τ (σ ) and τ (ρ) ) for larger m values and also depends on τ (ρ) , that is, the position of the peak in the spectrum relative to the analysed frequency range.

D I S C U S S I O N
The CCD-derived parameters m tot , τ mean and τ 50 increasingly diverge from the original CC model parameters with decreasing values of c and for τ (ρ) approaching the relaxation time limits in the decomposition.This result can be mostly attributed to the limited frequency range covered by the data, and thus the limited relaxation time range considered in the CCD.Outside the data frequency range the CCD response decreases faster than the original CC response (Figs 5a and c) (if the used kernel has a stronger frequency dispersion than the data), which causes a bias in the inferred mean and median relaxation times towards the centre of the RTD.However, a non-zero response outside the data range always remains due to the spectral width of the response of even a single kernel term, which can extend over more than one decade beyond the data frequency limits (increasing with decreasing c value, cf.Figs 5a and c).One decade can also be considered as the approximate theoretical resolution limit below which two peaks in a spectrum can no longer be distinguished (Florsch et al. 2014).Since the deviations of the CCD-derived relaxation times can easily exceed one decade, they must be considered significant and therefore should to be taken into account when analysing SIP data.
As notable in the RTD results in Fig. 5(b), in this example the CCDs for c = 1 (DD) and c = 0.5 (Warburg decomposition) yield similar chargeability weights at the data frequency edges, which For comparison the relationship between the corresponding CC relaxation times τ (σ ) and τ (ρ) according to eq. ( 4) is shown (dashed curves).
explains the larger spread of the response outside the data frequency range for the c = 0.5 case (Fig. 5a).However, with increasing original c value the contribution outside the data frequency range decreases, and thus the underestimation of m is also reduced (Figs 5c  and d).In fact for the trivial, but practically irrelevant, case of a decomposition kernel with identical frequency dispersion like the data (i.e.c = c), the derived parameters m tot , τ mean and τ 50 perfectly resemble the original CC parameters if the spectral peak lies well within the analysed data frequency range (cf.Figs 1-3).
The observed variation in the reconstruction quality of chargeabilities and relaxation times also influences the relationship between the mean relaxation times derived from the CCD in conductivity and resistivity formulations (Fig. 4).We found that τ (σ )  mean and τ (ρ)  mean agree only within one order of magnitude over the whole range of c values.Therefore, should τ (σ )  mean and τ (ρ) mean be compared or jointly used in a study, the discrepancy between both quantities should be taken into account.However, we again stress that errors can be much larger when CC and CCD-derived parameters are jointly interpreted and their systematic deviations ignored-no matter if conductivity or resistivity formulations are used.
We finally note that the CCD scheme depends on various settings which can have significant influence on the results, with the selected range of relaxation times having the strongest effect.As shown in previous work (Weigand & Kemna 2016), the relaxation time range should extend by at least one order of magnitude (better two) to both sides of the data frequency limits.The amount of regularization, here a smoothness constraint imposed on the RTD, has from our experience only marginal influence on the derived integral parameters m tot and τ mean .However, regularization at all is crucial to overcome the inherent ill-posedness of the problem (Florsch et al. 2014).

C O N C L U S I O N S
In this study we investigated the relationships of chargeability and relaxation time parameters recovered from a spectral decomposition of synthetic CC-type SIP data with the corresponding CC input parameters.As kernel function in the decomposition we considered CC models with different frequency dispersion, including the Debye Fidelity of spectral decomposition parameters 1419 model and the Warburg model.We found that total chargeability and mean relaxation time derived from the decomposition results can deviate significantly from their CC model counterparts, as well as among each other for different decomposition approaches.Total chargeability was underestimated by up to 80 per cent and mean relaxation time varied by up to three orders of magnitude compared to the CC input parameters.The deviations depend on the frequency dispersion of the analysed spectrum and the proximity of its peak to the frequency range limits considered in the decomposition.This behaviour is found for all considered kernels (Debye, Warburg, CC), except for the trivial, and practically irrelevant, case of a CC kernel with a frequency dispersion identical to the one of the data.Importantly, our results do not promote the use of a particular kernel function, which should be chosen based on the physics of the underlying polarization mechanism in the given application.
We conclude that a quantitative comparison of SIP parameters across different studies, or the adoption of parameter relationships from other studies, is only possible on the basis of a consistent, well documented spectral analysis procedure.This is especially important when comparing effective CC parameters with spectral parameters derived from decomposition results.Non-adherence to a consistent analysis procedure can result in severe misinterpretations in the growing field of applications, for example when linking SIP parameters to textural, hydraulic, or (bio)geochemical properties, or when transferring laboratory results to the field.

A C K N O W L E D G E M E N T S
Part of this work was conducted in the framework of the SFB TR32 'Patterns in soil-vegetation-atmosphere systems: monitoring, modelling and data assimilation' funded by the Deutsche Forschungsgemeinschaft (DFG).We thank the editor Jörg Renner and two anonymous referees for their very constructive reviews.

Figure 1 .
Figure 1.Relationship between original CC chargeability m and CCD-derived total chargeability m (ρ) tot using the indicated values of c for varying values of input CC time constant τ (ρ) and CC frequency exponent c (< c).A fixed value m = 0.5 was used in the CC model to generate the SIP data.The resistivity formulations of both the CC model and the DD scheme were used.

Figure 2 .
Figure 2. Relationship between original CC relaxation time τ (ρ) and the indicated DD-derived relaxation times for varying values of input CC parameters (τ (ρ) , c).A fixed value m = 0.5 was used for the input CC chargeability.Displayed is the log 10 -difference between τ (ρ) and the DD-derived relaxation times.The resistivity formulations of both the CC model and the DD scheme were used.

Figure 3 .
Figure 3. Relationship between original CC relaxation time τ (ρ) and CCD-derived relaxation time τ (ρ) mean for different values of c.The resistivity formulations of both the CC model and the DD scheme were used.

Figure 5 .
Figure 5. Debye versus Warburg decomposition of two CC model responses with different frequency dispersion.Decompositions are based on the shaded frequency range only (a,c).Original CC chargeability is m = 0.5 and CC relaxation time τ (ρ) = 0.159 s.(a) Responses of input CC model for c = 0.3 and CCD results for c = 1 (m (ρ) tot = 0.30) and c = 0.5 (m (ρ) tot = 0.35).(b) RTDs of the CCDs from panel (a).(c) Responses of input CC model for c = 0.5 and CCD results for c = 1 (m (ρ) tot = 0.41) and c = 0.5 (m (ρ) tot = 0.50).(d) RTDs of the CCDs from panel (c).The resistivity formulations of both the CC model and the DD scheme were used.