Debye decomposition of time-lapse spectral induced polarisation data

Spectral induced polarisation (SIP) measurements capture the low-frequency electrical properties of soils and rocks and provide a non-invasive means to access lithological, hydrogeological, and geochemical properties of the subsurface. The Debye decomposition (DD) approach is now increasingly being used to analyse SIP signatures in terms of relaxation time distributions due to its flexibility regarding the shape of the spectra. Imaging and time-lapse (monitoring) SIP measurements, capturing SIP variations in space and time, respectively, are now more and more conducted and lead to a drastic increase in the number of spectra considered, which prompts the need for robust and reliable DD tools to extract quantitative parameters from such data. We here present an implementation of the DD method for the analysis of a series of SIP data sets which are expected to only smoothly change in terms of spectral behaviour, such as encountered in many time-lapse applications where measurement geometry does not change. The routine is based on a non-linear least-squares inversion scheme with smoothness constraints on the spectral variation and in addition from one spectrum of the series to the next to deal with the inherent ill-posedness and non-uniqueness of the problem. By means of synthetic examples with typical SIP characteristics we elucidate the influence of the number and range of considered relaxation times on the inversion results. The source code of the presented routines is provided under an open source licence as a basis for further applications and developments. & 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Spectral induced polarisation (SIP) data consist of frequency dependent (typically 1 mHz-10 kHz) electrical impedances or admittances. These complex-valued data can either be represented in terms of real and imaginary parts, or as magnitude and phase values. Using the geometrical location of the electrodes, complex resistivities or complex conductivities can be inferred. These properties reflect electrical conduction and polarisation processes, for instance electrochemical polarisation at the interface between the pore fluid and minerals in soils and rocks or between different cells in living tissue, as of interest in geophysics and biophysics. SIP data are usually analysed and related to petro-, hydro-, or biogeophysical properties by means of empirical or mechanistic models (see, e.g., Slater, 2007;Kemna et al., 2012;Bücker and Hördt, 2013). While significant effort has been put into the development of mechanistic models over the last years (e.g., Revil and Florsch, 2010;Revil et al., 2012), most studies still address the formulation of parameter relationships using empirical models (e.g., Weller et al., 2013).
Two types of empirical model approaches can be distinguished: the first type describes SIP data using one or a few polarisation peaks (i.e., local maxima) in the absolute phase spectrum (or imaginary conductivity spectrum) with corresponding, distinct relaxation times. This includes the Debye model (Debye, 1929;Böttcher and Bordewijk, 1978) and the class of Cole-Cole-type models (Cole and Cole, 1941;Pelton et al., 1978, see Dias, 2000 for an overview). In general these models are overdetermined, containing far fewer model parameters than data points available. This, however, limits the flexibility regarding different shapes of the spectra and imposes strong data influence on fitting quality and robustness (Morgan and Lesmes, 1994).
The second approach describes the SIP response using a linear superposition of a large number of elementary Debye polarisation terms following a given distribution of relaxation times. The latter is commonly discretised in regular intervals. Thus, only the relative weights of the Debye contributions at discrete relaxation times are determined when fitting the observed polarisation response. This procedure of determining a relaxation time distribution (RTD) instead of a fixed number of relaxation times is referred to as Debye decomposition (DD). Contrary to the first type of models, the DD represents a strongly underdetermined inverse problem which requires some sort of regularisation to achieve a stable solution. A wide range of shapes of SIP spectra can be fitted using the DD approach, and more importantly, the approach does not require a certain number of polarisation peaks in the SIP data. However, representative (average or integral) parameters may be extracted from the RTD, also referred to as integrating parameters (Nordsiek and Weller, 2008) or characteristic integral parameters (Zisser et al., 2010b) in the literature.
The concept of the superposition of multiple polarisation terms goes back to Schweidler (1907), with the first proposition of a continuous RTD by Wagner (1913). One of the first published discrete decomposition procedures inverted the imaginary part of the complex permittivity (Uhlmann and Hakim, 1971), while later studies used the real part to fit permittivity spectra (Morgan and Lesmes, 1994;Lesmes and Morgan, 2001). It is noted here that because of the so-called Kramers-Kronig relationships (e.g., Macdonald, 1952) it is theoretically sufficient to only fit the real or imaginary part of the permittivity spectrum to infer the complete complex permittivity response (because both real and imaginary parts are related as a consequence of causality). Different implementations of the DD for time-domain data have been published (Tong et al., 2004;Tarasov and Titov, 2007), and likewise formulations for frequency-domain resistance and resistivity instead of permittivity or conductivity have been presented (Nordsiek and Weller, 2008;Zisser et al., 2010b;Florsch et al., 2012;Keery et al., 2012). Keery et al. (2012) investigated uncertainties of the inferred DD parameters using a Markov-chain Monte Carlo method. Recently, Florsch et al. (2014) proposed a decomposition approach based on relaxation models other than the Debye model, such as the Cole-Cole model or Warburg model. In an accompanying work, Revil et al. (2014) promote the use of the Warburg model, arguing that the elementary response in a rock is less dispersive than a Debye response. Following the same reasoning, already Tarasov et al. (2003) used a decomposition based on a Cole-Cole model in accordance with the response of a theoretical model for the polarisation of a single pore (Titov et al., 2002).
Based on numerous previous works on the linkage of SIP parameters with textural and hydraulic characteristics (for an overview see Kemna et al., 2012 and references therein), more recently also relationships of DD parameters to various petrophysical parameters have been established, such as to slag mass and grain size of slag-sand mixtures (Nordsiek and Weller, 2008), permeability of sandstones (e.g., Weller et al., 2010a), water saturation of sand-clay mixtures (Breede et al., 2012), specific surface area per unit pore volume of sandstones and sand-iron/clay mixtures (e.g., Weller et al., 2010b), and to changes in water salinity of sandstones (e.g., Weller et al., 2011). Attwa and Günther (2013) found correlations between DD relaxation time and hydraulic conductivity for samples taken from a Quaternary aquifer. One of the few published applications to field-imaging results investigated the relationship between DD-derived integral parameters and BTEX (benzene, toluene, ethylbenzene, and xylenes) contamination .
Despite the growing number of studies incorporating timelapse SIP data (e.g., Breede et al., 2012;Flores Orozco et al., 2013) for monitoring purposes, we are not aware of any DD approaches incorporating time-lapse data sets in a single inversion problem, for instance to apply regularisation with respect to temporal changes. In this work we present an open source implementation of a DD for complex resistance/resistivity SIP time-lapse (monitoring) data which is capable of imposing smoothness constraints with respect to spectral variations as well as changes over time. The inversion methodology is developed in detail, and integral parameters commonly extracted from DD results are briefly reviewed. Important aspects of the DD routine are examined, such as the choice of sampling (discretisation) of the relaxation time axis and the influence of regularisation strategies on the DD result. The source code (in the Python programming language) and corresponding documentation is available under an open source licence.

Methods
This section provides a detailed description of the DD methodology. Starting with the general formulation of the decomposition, the inversion algorithm is developed, including important implementation aspects, in particular the extension to invert timelapse SIP data sets. Fuoss and Kirkwood (1941) formulated the frequency response of the complex permittivity ε using a continuous relaxation time distribution, g τ ( ), which describes the Debye contribution at different relaxation times τ:

Debye decomposition
with ϵ 0 and ϵ ∞ being the asymptotic permittivities at low and high frequencies, respectively, τ the relaxation time, ω the angular frequency, and j 1 = − the imaginary unit. The Debye model follows from Eq. (1) for g 0 τ δ τ τ ( ) = ( − ), where δ denotes the Dirac delta function. In this case, τ 0 is inversely related to the frequency position (ω peak ) of the polarisation peak (in terms of imaginary part ϵ″) by 1/ peak The discrete case can now be formulated either by expressing g τ ( ) by a weighted sum of δ functions, or by discretising the continuous integral in Eq. (1). The discrete (and fixed) τ i values must sample at least the relaxation time range implicitly defined by the frequency range spanned by the data (i.e., according to the inverse correlation 1/ peak 0 τ ω = ) to facilitate an adequate approximation of the continuous case. SIP responses are now solely described by the weighting factors at the chosen relaxation times, representing a discrete RTD.
The analogon to the continuous formulation in Eq. (1) in terms of complex resistivity ρ ω( ) is given by with ρ 0 being the direct-current (DC) resistivity, N τ the number of relaxation times, m k the weighting factors, the so-called chargeabilities, corresponding to the relaxation times τ k . This formulation is based on the resistivity formulation of the Cole-Cole model by Pelton et al. (1978): ρ is the DC resistivity, m cc is the Cole-Cole chargeability, τ cc is the Cole-Cole time constant (which has the meaning of an effective relaxation time), and c cc is the frequency exponent, describing the strength of the frequency dependence (dispersion) (Fig. 1). Eq. (4) simplifies to the Debye model in terms of resistivity for c 1 cc = . It should be noted that the discrete form of Eqs. (1) and (3) are not fully equivalent to each other (Florsch et al., 2012;Tarasov and Titov, 2013), and care must be taken when comparing results obtained with the different formulations.
In this work, Eq. (4) is used to generate synthetic SIP data. Although in the literature the Debye decomposition has been mostly performed using the resistivity formulation (Eq. (3)), results in this work are presented in terms of complex conductivity, 1/ σ ω ρ ω( ) =^( ), which is commonly used as the basis for petrophysical interpretations (e.g., Revil and Florsch, 2010).

Integral parameters
From the RTD, different integral parameters can be derived, which provide information about the overall polarisation properties and relaxation time scales of the medium. We here consider the following parameters, which have been used in previous works: The total chargeability m m tot k N k 1 = ∑ = τ (Nordsiek and Weller, 2008) is the analogon to the Cole-Cole chargeability m cc . However, m tot only accounts for polarisations within the considered frequency range of the DD, whereas m cc also contains contributions from outside this range. Thus, not fully covered polarisation peaks imply an underestimation of the total chargeability m tot of the system.
The normalised total chargeability m m / n tot tot , 0 ρ = (Weller et al., 2010a) represents an extension of the normalised chargeability concept (e.g., Slater and Lesmes, 2002) to the total chargeability as derived from a DD. The normalised chargeability, m m/ n 0 ρ = , is commonly considered as a more appropriate measure of the strength of polarisation. In Appendix A we show for an individual Debye response according to Eq. (3) that the imaginary component of the complex conductivity, σ″, which represents a direct measure of polarisation according to established petrophysical models (e.g., Revil and Florsch, 2010), indeed directly scales with m n . For not too strong polarisations, i.e., m 1 ⪡ , one obtains Hence it is physically plausible to consider the normalised total chargeability derived from the DD in Eq.
(3) as a direct measure of the overall polarisation of the medium.
The cumulative relaxation time τ x denotes the relaxation time at which a certain percentage x of the total chargeability is reached (Nordsiek and Weller, 2008;Zisser et al., 2010b).
Commonly used is τ 50 , the median relaxation time of a given RTD.
The mean logarithmic relaxation time τ mean denotes the chargeability-weighted logarithmic mean value of the relaxation times (Nordsiek and Weller, 2008): As a generalisation of the maximum relaxation time τ max (e.g., Attwa and Günther, 2013), we here introduce the relaxation time l peak , τ , which refers to the lth local maximum of the RTD, with numbering beginning at low frequencies (i.e., high τ values). This parameter is useful if multiple dominant relaxation phenomena are to be analysed in a given SIP spectrum. These polarisation peaks result in multiple peaks in the RTD, which cannot be tracked using integrative parameters such as τ mean or τ 50 . The non-uniformity parameter U / 60 10 τ τ = τ is a measure of the width of the RTD (Nordsiek and Weller, 2008).
For a complete list of implemented integral parameters, we refer to the source code documentation.

Inverse approach
The inversion is implemented using a Tikhonov regularisation scheme (e.g., Zhdanov, 2002;Menke, 2012) which minimises an objective function, Ψ, comprising data misfit, Ψ y , and a model objective function, Ψ x , balanced by a real-valued regularisation where x is the model parameter vector, y is the data vector, and f x ( ) is the forward model response. We use individual data weighting factors which are stored in the diagonal matrix W . This corresponds to the assumption of uncorrelated data with Gaussian statistics, with the standard deviations of the data given by the inverse values of the diagonal entries of W . The regularisation term Ψ x is chosen to force smooth solutions via a roughness matrix R, which evaluates either the firstor the second-derivative variation in the components of x.
The inversion underlying the DD is formulated as a real-valued problem with real and (negative) imaginary components of the measured complex resistivity, obs i ρ ω ′ ( ) and representing the data at the measurement frequencies ω i . The model parameters are given by the log-transformed DC resistivity ρ 0 and the log-transformed chargeabilities m k at the considered relaxation times τ k . Our choice of using log-transformed chargeabilities in the parameterisation implies that relative changes in chargeability are equally weighted in the inversion, so that changes in the very low chargeability range still have importance in the inversion. Moreover, the log transform guarantees positive chargeability values. Ghorbani et al. (2007) suggest an alternative parameterisation of chargeability based on the logit function, which they argue is a mathematically more reasonable parameterisation also in the range towards the theoretical maximum value of chargeability (i.e., 1). However, since we are here not considering the total chargeability (like for instance represented by the Cole-Cole chargeability m cc ) but the chargeability contributions at the different relaxation times τ k , the values of m k are typically very small (i.e., 1 ⪡ ) and we therefore consider the log transform being used here adequate as well. Thus, the model parameter vector x, the data vector y (comprising one complex resistivity spectrum), and the forward model response f with N τ being the chosen number of relaxation times, N f the number of given frequencies, N x the number of model parameters in x N N 1 x ( = + ) τ , and N y the number of data points in y N N 2 The root-mean-square value of the misfit between data and model response for the imaginary parts (i.e., the second half of y and f x ( )), RMS Im , is used for assessing the quality of the inversion process: . 12 with the corresponding derivatives given in Eqs.
A cumulated sensitivity (or coverage), S Im k , , can be computed for each chargeability m k : This cumulated sensitivity value is a measure of the overall change of the imaginary parts of the predicted data in response to changes of an individual chargeability value m k , or in other words, how well the parameter m k is covered by these data. Correspondingly, large values of S Im k , indicate chargeabilities with a high influence on the response ρ″.

Model update
Minimisation of Ψ leads to an iterative Gauß-Newton like scheme, in which a model update x q Δ is computed for each iteration q using the normal equations: ]. The value of α is determined from a classical line search, explained further below, and the limitation 1 α ≤ prevents from overshooting (at the cost of potentially increasing the number of iterations until convergence is reached).

Starting model
The starting model for each spectrum is determined by testing a series of homogeneous RTDs for chargeability values between 10 À 12 and 1. The model yielding the minimal value of RMS Im is then used as the starting model in the inversion.

Stopping criteria
The inversion is stopped when one or more of the following criteria are met: The RMS Im does not change significantly (smaller than predefined percentage) between subsequent iterations.
The RMS Im decrease lies below a certain threshold value. An exception is the first iteration, where an RMS Im increase of up to 100% is allowed to account for a possible tuning-in phase based on the starting model.
The update leads to a numerical error (e.g. overflow error).

Step-length selection
At each iteration an optimal step length parameter α (Eq. (16)) is selected using a line search approach (e.g., Kemna, 2000;Günther, 2004): a parabola is fitted through three RMS Im values corresponding to the α values 0, 0.5, and 1, and the minimum of the parabola determines the step length α min for the model update. In case of a minimum above 1 the step length is set to 1. For 0 min α ≤ , the model update does not lead to an improvement of RMS Im and the inversion is stopped.

Data weighting
The diagonal entries of W compensate for different ranges of ρ ω ′( ) and ρ ω ″( ): This weighting scheme ensures that both real and imaginary parts are fitted to similar misfit levels irrespective of their values.

Regularisation
Smoothness constraints are applied to the chargeability values using firstand second-order finite-difference operators (e.g., Menke, 2012;Aster et al., 2013). The first-order smoothing operator is given by the

Selection of τ values
Two parameters determine the discrete relaxation times used for the decomposition: The number of relaxation times per frequency decade, N d , and the min/max limits for the relaxation times. The N d parameter must provide enough degrees of freedom to fit the data (Uhlmann and Hakim, 1971). Minimal and maximal τ values can be inferred from the frequency limits of the data using the inverse correlation of τ and ω (Uhlmann and Hakim, 1971).
These limits can be extended by multiplicative constants, defaulting to an extension of one frequency decade to either side of the data frequency limits. Integral parameters are only computed for the τ range defined by the data frequency range, irrespective of the actual range set for the decomposition.

Selection of λ freq values
Three methods of selecting the regularisation parameter are implemented: a fixed value can be used for all iterations, the Lcurve (e.g., Hansen, 1990;Florsch et al., 2014) can be used to infer an optimal regularisation parameter (Zisser et al., 2010a;Florsch et al., 2012), and a line search (simplified from the one used in Kemna, 2000) is implemented to automatically determine an optimal λ freq value at each iteration. In the latter approach, the λ freq value of the previous (qth) iteration is used as a starting point and multiple new values in the range between 0.1 freq λ and 10 freq 4 λ are sampled. The value yielding the smallest RMS Im value is then used to actually compute the model update. For the first iteration, the λ freq value is set to the number of model parameters. The inversion results presented in the following were computed using this last approach of selecting λ freq .

Extension to time-lapse inversion
Time-lapse data can be accounted for in the inversion scheme by extending the vectors and matrices to include all time steps. Based on the formulation for a single spectrum ρ ω( ) the new vectors sequentially contain the quantities for the different time steps, and the new matrices become block matrices with each block along the diagonal corresponding to a time step. In the following, the subscript d denotes the time index and N t the total number of time steps.
The The corresponding Jacobian matrix can be assembled based on Eq. (13) for each time step: The RMS Im all value incorporates the misfits of the imaginary parts for all time steps: ,  (Eq. (19)). The smoothness constraints along the time axis require some further consideration. The smoothing operator has to be applied to each parameter of x separately along the time axis, and thus correspondingly resized versions of Eqs. (18)  The scheme described above simplifies the implementation of the regularisation matrices, as only a base version (with variable size) has to be created for each regularisation type, which then can be applied to both frequency and time domains using the projection. It is noted that so far the smoothing operators do not take the distance between adjacent parameters in the respective direction (frequency or time) into account. However, in particular for timelapse measurements it is not uncommon to encounter irregular time intervals t Δ . These time intervals can be used to weight the regularisation matrices by means of a diagonal time-weighting matrix C t : with subscript tw indicating the time-weighted variant of the regularisation matrix. The time weighting can be disabled by setting the matrix C t to the unity matrix. In the current implementation, time weighting can only be applied in conjunction with first-order smoothing. However, in the future we plan to also extend it to second-order smoothing based on the implementation of the second-order derivative operator also for non-regularly sampled points.

Results and discussion
In this section the optimal selection of the number and range of τ values is examined, followed by an investigation of the effect of the regularisation, both in frequency and time domain.

Relaxation times
The choice of relaxation times is crucial for the decomposition process. Using an exemplary bimodal SIP response (Fig. 2a), recovered integral parameters were investigated for various N d values (in the literature, values for N d range from 4 to 166 (Uhlmann and Hakim, 1971;Morgan and Lesmes, 1994;Nordsiek and Weller, 2008;Zisser et al., 2010b).  The response of even a single Debye term spreads significantly over a range of more than one frequency decade to both sides of its peak frequency (e.g., Florsch et al., 2014). Hence, τ values outside the data frequency range can influence the inverted RTD inside the data frequency range. To illustrate the effects of too narrow τ ranges, synthetic SIP data were generated in a frequency range between 10 À 3 Hz and 10 Hz 3 using a two-term Cole-Cole model (Fig. 3a), and subsequently the DD was applied using three different τ ranges: (1) the range was determined from the data frequencies, (2) the data frequency range was symmetrically extended by one frequency decade, and (3) the data frequency range was symmetrically extended by two frequency decades. In all cases, N d was set to 20, and the data could be adequately fitted. However, the obtained RTDs vary considerably (Fig. 3b-d). Without any extension of the τ range (Fig. 3b), two peaks can be observed in the RTD (Fig. 3c), of which only one can be related to a polarisation peak in the SIP data (compare detected peaks, i.e., solid and dashed lines in Fig. 3a and b). However, when the τ range is increased (Fig. 3c and d), only the polarisation peak indeed present in the data is recovered. The high-frequency peak (low relaxation times) in the RTD (Fig. 3b) corresponds to a large increase in S Im (Fig. 3e), indicating an over-proportional influence of the small τ values in the case of the insufficiently broad τ ranges. Note that, irrespective of the τ range used, large values of S Im can be observed for the not fully covered high-frequency peak, suggesting a strong influence of the corresponding chargeabilities if polarisation peaks are not fully covered by the data.
Based on the results (Figs. 2 and 3) we recommend to use a value of at least N d ¼20 for the DD, and to extend the τ range at least by one frequency decade to each side of the data frequency limits.

Regularisation strategies
The effect of too weak or too strong (first-order) frequency regularisation is demonstrated using synthetic SIP data contaminated by normally distributed noise (this type of noise is in the following also referred to as measurement noise). Three DDs were performed using regularisation parameter values 1, 10 3 λ = , and 10 5 . The corresponding fitting results are shown in Fig. 4a, c, and e, with RTDs plotted in Fig. 4b, d, and f. For a λ value of 1, a relatively rough RTD is obtained, nonetheless adequately fitting the data. Due to the roughness of the RTD no clear peak relaxation time can be identified. The DDs performed with λ values of 10 3 and 10 5 show smooth RTD curves, each exhibiting a clear peak. However, the larger regularisation parameter leads to a flat polarisation response which does not fit the SIP data (Fig. 4e). Thus, the used regularisation parameters should be verified to avoid situations as presented in Fig. 4b and f.
Considering the analysis of time-lapse data, no inherent constraints (or preferences) regarding the time evolution of the DD parameters can be inferred from the method itself. All such constraints must be based on some sort of a-priori information (e.g., on the type of polarisation or the type of process being investigated with SIP measurements). Therefore, a variety of regularisation operators can be used for time regularisation. However, it should be kept in mind that different constraints reflect different a priori information, which can have considerable influence on the Fig. 3. DD results using different τ ranges. (a) Synthetic SIP data (dots) and fitting results (solid curves) (grey: σ′, black: σ″). (b,e): RTD and coverage S Im using τ limits determined by data limits. c,f): RTD and coverage S Im using τ limits extended by one frequency decade (shaded area) relative to the data limits. (d,g) RTD and coverage S Im using τ limits increased by two frequency decades (shaded area) relative to the data limits. All fits were performed using a frequency regularisation parameter 50 λ = and N d ¼ 20. Vertical lines mark peaks in the inverted RTD (see also discussion in the text). resulting RTD.
The influence of time regularisation operators on the decomposition process is illustrated on a synthetic time-lapse SIP data set which comprises spectra at 20 (irregular) time steps created using a single-term Cole-Cole model response (Fig. 5). The chargeability values of the SIP spectra decrease linearly with time, imitating a corresponding time-dependent polarisation process (Fig. 5a, increasing opacity indicates increasing time). A small, normally distributed measurement noise component with a standard deviation of 0.5 mrad was added to each spectrum. Additionally, a time dependent, normally distributed noise component with a standard deviation of 5% was added to the Cole-Cole chargeability values to simulate variations in the underlying polarisation process (in the following also referred to as process noise). This second noise component does not distort the SIP spectra, however, the spectra are vertically shifted. This is the primary type of noise that can be smoothened out using a timeregularisation approach. Measurement noise would have to be very large in order to influence integral parameters such as m tot . The m n tot , values inferred from the inverted RTDs for different time-regularisation approaches in the DD are plotted in Fig. 5b. For consistency, the DD without any time regularisation was performed with the same frequency regularisation parameter as the DD with time regularisation. The results show that time regularisation successfully smoothens the noise-induced roughness in the m n tot , evolution compared to the DD results without time regularisation. Due to the non-regular spacing of the time-steps (as often found in long-term experiments, e.g. due to equipment failures or miscellaneous delays), the time-weighted regularisation scheme creates smoother results in this example.
Time regularisation operates on the RTD, i.e., the chargeabilities at the discrete relaxation times. Smoothing of the chargeabilities thus directly influences the integral chargeability parameters m tot and m n tot , , as shown above. Integral relaxation time parameters are only influenced implicitly due to the deformation of the RTD by the time regularisation. Nonetheless, time regularisation can improve the reconstruction quality of integral relaxation time parameters for noisy time-lapse SIP data. Here, measurement and process noise components are treated separately. Fig. 6a presents 20 single-term Cole-Cole model responses whose relaxation time changes linearly with time on a logarithmic scale (the simulated time steps, again, are not regularly spaced). A normally distributed process noise level Fig. 4. DD results for different values of the frequency regularisation parameter λ freq . a,c,e) Synthetic SIP data (dots) and fitting results (solid curves) for λ freq values of 1, 10 3 , and 10 5 , respectively (grey: σ′, black: σ″). Uncorrelated, normally distributed noise with standard deviations of 1.5 m Ω and 2 mrad for magnitude and phase of the complex resistivity, respectively, was added to the SIP data prior to DD. (b,d,f)    of 20 % was added to the log-values of the relaxation times before computing the responses. This large noise level is justified by the large possible dynamics that can be observed for relaxation times. Additionally, a small measurement noise level of 0.5 mrad was added to the generated spectra. The time evolution of the integral parameters peak 1, τ and τ mean , inferred from the obtained RTDs, is presented in Fig. 6b and d, respectively. Time regularisation greatly improves the quality of the τ mean results, while the peak 1, τ results are only slightly improved. This can be explained by the nature of the data noise: process noise components in τ directly translate to changes in the peak position of the RTD, without affecting the smoothness of the spectrum and the RTD. Time regularisation strategies thus have to move the RTD peak in order to change peak 1, τ results, but only have to change the shape of the RTD in order to influence τ mean . This is illustrated in Fig. 6c and e. Time regularisation changes the model response only slightly, but deforms the RTD enough to influence τ mean (Fig. 6c, black curve).
However, peak 1, τ results can be heavily improved by time regularisation schemes, if the data are primarily influenced by measurement noise. This effect is illustrated in Fig. 7, which shows a simulation similar to the previous one (presented in Fig. 6). However, no process noise component was added to the Cole-Cole relaxation time, and the standard deviation of the measurement noise was increased to 2 mrad. As can be observed in Fig. 7b, peak 1, τ is heavily influenced by these noise components, while τ mean is only slightly affected. This result can be explained by the resulting RTDs (Fig. 7e), where the noise produces misdirecting peaks which get smoothed out by the time regularisation.
Note, again, that peak 1, τ is usually only interesting if multiple dominating polarisation responses are present in the SIP spectra, and only a subset is analysed. In other cases a more robust relaxation time parameter such as τ mean should be preferred.
The presented results  illustrate that care must be taken with respect to both adequately balancing regularisation versus data misfit and adequately balancing time versus frequency regularisation. The former issue can be addressed by, for example, performing a line search at each iteration of the DD to find the regularisation parameter value which locally minimises the data misfit, or by using the L-curve criterion. For the present inverse problem, however, it was found that the latter approach may yield varying results. To ensure a good balance between frequency and time regularisation, the involved regularisation parameters should be chosen such that the norms of the different regularisation terms in Eq. (26), i.e. R x R x , f all all 0 ρ , and R x m all , yield values of similar order of magnitude, so that the dominance of a single term over the others is avoided.

Conclusions
The presented Debye decomposition procedure is a robust and versatile tool to describe and analyse time-lapse SIP data. The implementation of the procedure available under an open source licence facilitates further development and an easy adaptation to specific data types and fields of application. Time-regularisation strategies yield improved results for a better interpretation of monitoring measurements, especially when peak relaxation times are considered. The number of relaxation times per frequency decade should be larger than 20, and the range of relaxation times should be extended by at least one frequency decade beyond the limits given by the data frequencies. In addition, as with all regularised inversion problems, the choice of regularisation type and strength can significantly influence the results of the Debye decomposition. Future developments of the DD implementation will include the incorporation of specific polarisation terms, for example to represent higher-frequency electromagnetic coupling responses, as well as the inclusion of spatial regularisation for an improved analysis of SIP imaging results.
The code described in this study is maintained and developed at https://github.com/m-weigand/Debye_Decomposition_Tools.