Effect of current pulse duration in recovering quantitative induced polarization models from time-domain full-response and integral chargeability data

The signal level and shape of induced polarization responses are signiﬁcantly affected by the current pulse duration and waveform. If not accounted for, this data dependency on the current will propagate trough the inversion to results rendering unquantiﬁable subsurface models. While this problem has been addressed in full-response induced polarization modelling, questions remain as to how to accurately retrieve quantitative induced polarization inversion models from the types of apparent integral chargeability data often used in data interpretation. Although several methodologies have been proposed for handling and inverting apparent resistivity and integral chargeability, these cannot compensate for the data dependency on the current waveform and pulse duration. This paper presents a novel inversion method for such data. The method considers current waveform and receiver transfer functions for retrieving quantitative IP models unbiased by transmitter waveform. The method uses the constant phase angle model, expressed in terms of the medium resistivity and phase. Speciﬁcally, four ﬁeld data sets for the same proﬁle but with different 100 per cent duty cycle pulse durations (4, 2, 1 and 0.5 s) serve as examples of data sets giving models dependant on current waveform when inverted with standard approaches. The novel inversion method presented here gives quantiﬁ-able models independent on the current waveform and pulse duration. These results resemble models retrieved with existing, full-response induced polarization inversions. The results still


I N T RO D U C T I O N
Previous studies have pointed out various consequences of varying current pulse duration in time domain induced polarization (IP) measurements. The current waveform, especially the pulse duration, influence the IP data (Van Voorhis et al. 1973;Fiandaca et al. 2012Fiandaca et al. , 2013Olsson et al. 2015b;Mao et al. 2016) and the signal-to-noise ratio increases with increasing pulse duration Olsson et al. 2015b). Studies of inversions that do not consider the waveform have found that they can only provide qualitative IP information for the subsurface (Fiandaca et al. 2012Olsson et al. 2015b;Mao et al. 2016). Parameter studies (Lajaunie et al. 2016) and Markov chain Monte Carlo methods (Madsen et al. 2017) have investigated effects on the spectral information content and its dependence on pulse duration. However, the question of how to retrieve quantitative IP inversion models independent of current pulse duration from time-domain apparent integral chargeability data has received little attention.
Several inversion methods exist for resistivity and induced polarization data from time domain measurements and some of these handle apparent integral chargeability data. Briefly, the four main methodologies are: (1) inversion for chargeability with two DC forward calculations, one for DC and the other for DC modified with the chargeability (Oldenburg & Li 1994), (2) complex resistivity inversion for time domain data transformed into frequency domain (Kemna et al. 2000), (3) complex resistivity inversion with one separate DC forward calculation modified with chargeability for each IP gate (Hönig & Tezkan 2007) and (4) complex resistivity full-response inversion with modelling of transmitter and receiver waveform (Fiandaca et al. 2012. Quantitative inversions of full-response IP data (with multiple gates per IP response) based on Cole-Cole parametrization have been previously retrieved using the fourth methodology listed above (Fiandaca et al. 2012), but such highly parametrized models (four model parameters) are not suitable for single IP datum. This is the case for apparent integral chargeability data. However, the methodology has also been successfully applied with constant phase angle (CPA) parametrization (Johansson et al. 2015;Olsson et al. 2015b), an approach that uses only two model parameters. Approaches with fewer model parameters could be suitable for apparent integral chargeability data. Manual processing for outliers in full-response IP data is at present much more time consuming than for integral chargeability data. In many cases, this former approach is used for rapid post-processing and inversion procedures, namely for getting preliminary results in the field. The fact that full-response IP data at present does not qualify as industry standard data and the long-term ubiquity of apparent integral chargeability data require development of a methodology that can handle apparent integral chargeability data in a manner unbiased by transmitter waveform and pulse duration.
This paper presents a novel methodology for inverting apparent resistivity and apparent integral chargeability IP data based on integral CPA (iCPA) modelling. This latter inversion approach reduces the effect of current waveform and pulse duration on inversion models. The inversion methodology is compared to two of the existing inversion methodologies (2 and 4 above) using data sets with varying pulse duration acquired along a 2-D profile. Specifically, four field data sets for the same profile but with different 100 per cent duty cycle pulse durations (0.5, 1, 2 and 4 s), spanning almost a decade in duration, serve as examples of data that generate unquantifiable inversion models when inverted using the standard methodology. We show that the iCPA methodology clearly reduces effects related to pulse duration on inversion models. However, we also show that parametrization uncertainty remains for the iCPA modelling due to the assumption that the data can be accurately described with a CPA parametrization. This remaining parametrization uncertainty is further discussed based on full IP response CPA and maximum phase angle (MPA, Fiandaca et al. 2018) modelling for the same data sets.

Induced polarization field data
A time domain resistivity and induced polarization field survey was conducted at the construction site for the European Spallation Source (ESS), Lund, Sweden in 2014. The local geology at the site consists of clayey till overlying Silurian shale with dolerite dyke intrusions. The survey profile was acquired with 64 acid-grade stainless steel electrodes spaced at 2.5 m to give a profile length of 157.5 m. The profile was approximately centred over and perpendicular to a known dolerite dyke serving as a chargeable IP anomaly. Two separated multiconductor cables, one for transmitting current and the other for receiving potentials, served to reduce wire-to-wire coupling and improve data quality (Dahlin & Leroux 2012). To enable dense measurements, the cables used for current transmission and potential measurements were swapped after measuring the electrode combinations for one cable arrangement. An ABEM Terrameter LS instrument was used for transmitting current, measuring potentials and recoding full waveforms of these data at 1 kHz sampling rate. The electrode contact resistance was estimated for all electrodes with the Focus One protocol (Ingeman-Nielsen et al. 2016) resulting in a mean contact resistance of 210 and a standard deviation of 90 . Four field data sets of 1015 measurements each were acquired on the same profile using a 100 per cent duty cycle current injection (Olsson et al. 2015a) of four stacks with different pulse durations. These durations were set at 0.5 s, 1 s, 2 s and 4 s. The recorded full waveform data were subsequently processed for harmonic denoising, despiking and background drift removal following Olsson et al. (2016). The IP responses were gated with approximately logincreasing IP-gates with the same temporal distribution but with more gates for the longer on-time acquisitions. The processed data sets were manually inspected for obvious outliers in the software Aarhus Workbench (version 5.6.3.0) which allows the analyst to individually filter the gates of the full IP response. The integral chargeability data used for the standard inversion were processed in Res2Dinvx64 (version 4.06.07) for obvious outliers in integral chargeability data measured from gate 5 to 19. These data approximately corresponded to IP response observations between 5 and 500 ms. Fig. 1 shows the resulting fully processed data as pseudo-sections of apparent resistivity (top) and apparent integral chargeability for the four different pulse durations. All data sets show smooth variations after processing, which indicate that the data are generally of good quality. The figure shows one apparent resistivity pseudosection in absolute magnitude (the 4.0 s data set), while the pseudosections of the other pulse durations are instead shown as relative differences compared to the 4.0 s apparent resistivity, in order to highlight the data differences. As expected based on previous research (Mao et al. 2016), there is a general trend of decreasing magnitude in apparent resistivity data with decreasing pulse duration and the effect is larger for datapoints with higher apparent chargeability. The apparent integral chargeability data sets vary significantly with varying pulse durations. Magnitudes increase with increasing pulse duration but variations within the same data set are similar with higher magnitudes at greater depths. This effect arises from longer pulse duration which contribute more energy into the subsurface thereby charging a wider spectrum of polarization processes (Mao et al. 2016). Fig. 2(left) shows individual IP responses for the same quadrupole from the four different data sets corresponding to an apparent focus position of x = 65 m and z = 16 m (magenta cross in Fig. 1), its corresponding full waveform for the 4.0 s pulse duration (mid) and the full waveform for a shallower measurement with modest IP effects seen in the data (of x = 124 m and z = 5 m, blue cross in Fig. 1). As expected, longer pulse durations increase the measured chargeability (in mV V -1 ). Furthermore, the IP response for all gates show this effect. This pulse duration dependency clearly demonstrates that the inversion process must consider the current waveform and pulse duration in order to retrieve quantitative IP models from time domain measurements. Furthermore, the comparison of full-waveform data for quadrupoles with high-and low-chargeability values highlights the IP effect on the DC data, with the high-chargeability full-waveform data still increasing the potential values at the end of the pulses.

Standard inversion of integral chargeability data
Standard inversion for the second methodology has been carried out in Res2Dinvx64 (version 4.06.07). The inversion used default inversion settings with the exceptions of a lowered stopping criteria limited to a relative RMS error change of 1 per cent, an increase  in maximum number of iterations to 15, the use of extended model (rectangular) and an increase in model discretization depth to approximately 54 m, with 24 discretized layers and increasing thickness factor of approximately 1.117. All inversion models converged after 13-14 iterations with an RMS error of 0.62-0.64 per cent for resistivity data and 0.46-0.96 per cent for IP data. Fig. 3 shows the resulting inversion models for the IP data sets shown in Fig. 1. As expected, the inversion of different data sets gives highly similar resistivity inversion models that exhibit a low resistive top layer corresponding to the clayey till and an underlying unit of higher resistivity corresponding to the Silurian shale. The IP models show a chargeable anomaly at the position of the dolerite dyke but give different magnitudes for different data sets, with higher magnitude for data sets with longer pulse durations. The IP magnitude variation arises from the apparent integral chargeability dependency on the current waveform. This shows that the standard inversion methodology cannot quantitatively recover IP material properties independent of current waveform.
For obvious reasons, geophysical methods need to retrieve consistent inversion models for the same survey profile independent of acquisition settings in order to accurately estimate material properties of the subsurface. Previous approaches have attempted to include the full IP response in the inversion and model the current waveform based on Cole-Cole (Fiandaca et al. 2012 or CPA (Johansson et al. 2015) parametrizations.

Inversion of full response IP data
Full IP response inversion using the fourth methodology (complex resistivity full-response inversion with modelling of transmitter and receiver waveform) has been carried out in Aarhusinv . The input data for this inversion consists of the apparent resistivities, ρ a , and all the individual IP gates after post-processing, M i , giving data space d obs = {ρ a , M i } for i = 1:N i where N i represents the number of gates.
The model space of the inversion is defined through a parametrization of the complex resistivity spectrum, assigning the parametrization values to the cells of the modelling 2-D mesh.
The objective function of the full-response complex resistivity inversion is defined as: where δd is the data misfit, C obs is the covariance matrix of the observed data, δr is the model roughness and C R is the covariance on the roughness constraints. N d and N R are the numbers of data parameters and roughness constrains, respectively. The objective function is minimized through a gradient-based 2-D inversion scheme that is described in detail by Fiandaca et al. (2013) and Auken et al. (2015). For the inversion, a noise model was used with 1 per cent relative noise in DC and 10 per cent relative noise in IP together with a voltage dependant noise threshold of 0.05 mV for a nominal integration time of 0.01 s and one stack (Olsson et al. 2015a). Vertical and horizontal model roughness constraints of 1.3 and 1.2, respectively (expressed as standard deviation factor or STDF) were applied to the full model space for the inversion. The model space was discretized into 22 vertical layers with log-increasing layer thickness down to a maximum depth of 52 m and parsed horizontally into 2.5 m columns.
The parametrizations used in this study are based on two complex resistivity models, that is the Drake model (Van Voorhis et al. 1973) and the resistivity Cole-Cole model (Pelton et al. 1978), as described in eq. (1) and (2), respectively: In eq. (2) K is a positive constant, ω l represents the low-frequency pole that allows for the definition of the DC resistivity ρ 0 = K(ω l ) −b and b is proportional to the phase of the complex conductivity ϕ = π 2 b for ω > ω l . In eq. (3) ρ 0 represents the DC resistivity, m 0 is the intrinsic chargeability, τ ρ is the time constant (defined as the inverse of the angular frequency at which the imaginary part of the complex resistivity has a maximum) and C is the frequency exponent.
In this study, instead of using the parameters of eq. (2) in the inversions, the constant phase angle (CPA) parametrization is used, in terms of the parameters {ρ 0 , ϕ} with fixed ω l = 10 −6 rad/s (Johansson et al. 2015). Furthermore, the parameters {ρ 0 , ϕ max , τ ϕ , c} of the maximum phase angle (MPA) model ) are used instead of the classic Cole-Cole ones {ρ 0 , m 0 , τ ρ , c}. In the MPA model m 0 is replaced by the maximum phase ϕ max of the complex conductivity, and the time constant τ ρ by τ ϕ , that is the inverse of the angular frequency at which ϕ max is reached. The MPA re-parametrization effectively reduces parameter correlations and provides better resolved inversion models . Furthermore, for low values of the frequency exponent C the maximum phase angle ϕ max of the MPA model tends to the constant phase ϕ of the CPA model, making MPA and CPA models equivalent for low C values.
For the full-response IP inversion, this study used the CPA parametrization as a starting point for inversion. This starting point is motivated by similar shapes for most of the data as plotted in log-log space (with different magnitude), which indicate that the relatively simple CPA model in this case could suffice to explain a majority of the data (Lajaunie et al. 2016). Fig. 4 shows CPA inversion models for resistivity (top row) and CPA phase angle (second from top row) together with the estimated depth-of-investigation (DOI, Fiandaca et al. 2015) for STDF thresholds of 2 and 5 for the shallower and deeper DOI. Fig. 4 also shows the apparent resistivity (third from top) and misfit pseudo-sections of apparent integral chargeability (fourth from top). Focus positions marked in the pseudo-sections indicate two measured IP responses along with their corresponding forward responses (bottom), corresponding to the full waveform data shown in (Fig. 2). All inversion models converged after 9-13 iterations with a final chi misfit of between 0.64 and 0.76 for DC and between 0.90 and 1.1 for IP. Both the resistivity and CPA phase angle inversion models are similar for all four pulse durations. Notably, the inversion models show a highly resistive anomaly in the position of the dolerite dyke whereas the resistivity models retrieved with the standard methodology are homogeneous below about 70 masl. This arises from the DC potential correction (Fiandaca et al. 2012 which assigns the dyke greater resistivity than that estimated by the standard inversion when high polarization effects are present, due to the fact that the potential does not reach the DC asymptotic value at the end of the pulses, as evident in Fig. 2. It could also arise from contrasting resistivity equivalences in the respective software packages since inclusion  of full response IP data in the inversion is known to reduce problems associated with resistivity equivalences . The phase angle models show a highly chargeable anomaly (above ∼50 mrad) occurring in the upper left part of the main chargeable anomaly. The anomaly occupies the position of the dyke, but the phase anomaly magnitude increases with longer pulse durations. In general, some IP measurements persist with relatively high IP misfit, especially in upper regions of the profile between 10-60 m and 110-150 m. An example of the data (blue) and forward response (overlying black) in the latter coordinate interval is shown in the figure (bottom) and exemplifies the elevated misfits. Note the large difference between the 4.0 s pulse duration forward response and associated data relative to the deeper response (magenta). The high and systematic misfit of the blue response indicates that the CPA model is not adequate under these circumstances. Alternatively, it may indicate that the IP datum is erroneous due to poor background drift removal, since upper regions have low polarizability giving low IP signal levels. The overall IP misfit also decreases with decreasing pulse duration, as expected for subsurface images from a Cole-Cole based model but inverted with a CPA model. This occurs due to the shorter IP response time range (i.e. pulse duration) making the Cole-Cole IP responses more CPA-like (Lajaunie et al. 2016). The blue response gives rise to same effect of lower misfit for shorter pulse durations. Inversions such as the MPA, described here, use additional parametrization for increasing the forward response flexibility. Fig. 5 shows MPA inversion models for resistivity (top) and the MPA Cole-Cole model parameters (second to fourth) together with the estimated DOI  as STDF. The models are shown with ρ 0 and ϕ max STDF values of 2 and 5, 10 and 25 for τ ϕ and 1.5 and 3 for c. The use of different STDF thresholds for the different parameters attempts to reflect the variability in each parameter magnitude. Fig. 5 also shows the apparent resistivity (fifth) and apparent integral chargeability (sixth) misfit pseudo-sections and examples of two observed and forward calculated IP responses (bottom) with their focus positions marked in the pseudo-sections (same position as Fig. 4). All inversion models converged after 9-11 iterations with a final misfit of between 0.63 and 0.73 for DC and between 0.89 and 1.0 for IP. Further parametrization of the model using a MPA model improved the IP misfit slightly for high misfit areas, as evidenced by the high misfit response (blue IP response, Figs 4 and 5). However, in general the MPA model offered only marginal improvement on the overall IP data misfit giving a maximum mean misfit difference of approximately 0.1. The longest pulse duration, which exhibits the highest spectral information content, gives the largest difference. This supports the interpretation that the CPA model does not suffice in explaining all the IP data. The MPA model provides better results for the IP responses in the upper part of the pseudo-section, which are dominated by the contribution from the clayey till. Areas with higher misfits remain however, especially in upper areas with lower IP signal levels. Misfit is often systematic and within error envelopes. These results indicate that the MPA model, despite being an improvement on the CPA model, still fails to precisely capture the spectral content but this can also be an indication of systematic noise in the data, perhaps due to insufficient removal of background potential drift.
In general, this methodology, which is based on full-response IP inversion and modelling of transmitter and receiver transfer functions, eliminates the effect of pulse duration seen for the standard methodology (Fig. 3) and gives similar phase magnitude IP models independent of pulse duration.

Novel inversion of integral chargeability data
A novel method for inverting apparent integral chargeability data that utilizes the CPA parametrization for forward modelling of the waveform dependant apparent integral chargeability has been developed and implemented in the software used for the full response inversions. Similarly to the standard inversion of integral chargeability, the data space for this inversion scheme, d obs = {ρ a , M int }, consists of apparent resistivity ρ a and apparent integral chargeability: where M i and t gate,i are the chargeability value and time width of the ith gate and t tot = N i i=1 t gate,i is the total time summed from the first to the last non-filtered gate. This means that in the integral chargeability CPA (iCPA) approach, integration interval can vary for different quadrupoles within the same data set depending on data quality. Thus, it is possible to do adaptive rejection of individual IP gates and to keep the noise free subparts of noisy responses. In the iCPA implementation, the STD for each IP response gate M i , STD i , is propagated to the M int datum as uncorrelated error: The forward routine is based upon existing procedures for frequency to time transform and modelling of received potential stacking and gating for decoupling IP inversion models from current waveforms (Fiandaca et al. 2012), similar to methods used in standard full response inversions. However, the forward response computation is modified for retrieving the apparent integral IP. The modifications involve splining of the time domain response, specifically for example, by calculating for ten points per decade. The analytical integral for the spline is calculated from the first to the last active IP gate in the forward data. This method for handling integral chargeability data has been implemented for integral CPA (iCPA) modelling, for both the 50 per cent and the 100 per cent duty cycle current waveforms. Fig. 6 shows the resistivity (first) and iCPA phase angle (second) inversion models together with the estimated DOI calculated for STDF thresholds of 2 and 5 as well as apparent resistivity (third) and apparent integral chargeability (forth) misfit pseudo-sections. The data sets correspond to those used for the standard methodology and full-response inversions but with the IP responses summed between the first and last non-filtered gate. The iCPA inversions were carried out in the same manner as the CPA inversions. They specifically used the same settings and discretization, but a data space set according to eq. (4). The iCPA inversion also used the same noise model as that used in the full response CPA inversion for DC and IP data but propagated the individual IP gate errors to a single error for the apparent integral chargeability datum according to eq. (5). All inversion models converged after 11-14 iterations with a final chi misfit of 0.59-0.64 for DC and 2.2-3.3 for IP. The resistivity inversion models for the different pulse durations are similar to each other and similar to the CPA and MPA inversions. The DOI for the CPA resistivity and phase angle inversions give generally deeper results compared to the iCPA inversion models, reflecting the lower overall chi misfit.
The misfit pseudo-sections show evenly distributed misfits except for higher misfits in shallow areas and towards the edges of the profile. The overall IP misfit is generally elevated relative to the full response inversions. The IP misfit reaches mean values of 2.2 and maximum value of 3.3 indicating that the CPA model is not capable of explaining all IP responses and that the misfit is systematic. However, for the sake of direct misfit comparison between the iCPA inversion model and full-response IP inversion models, the integral chargeability misfits of the full response inversions are estimated by summing the data and forward responses according to eq. (5)

ST D i t gate,i .
(6) Table 1 shows a summary of the propagated integral apparent chargeability misfits for the CPA and MPA inversion models as well as the iCPA inversion model misfits. The uncorrelated apparent integral chargeability misfits for the CPA and MPA inversions are in fact within the same range as the iCPA inversion model misfits when each data set is considered independently. Additionally, if the STDs are propagated as correlated error, all misfits fall below one for CPA and MPA results. This indicates that the iCPA, CPA and MPA inversions give equally valid results (while the MPA gives slightly better results). Remaining elevated misfits arise from correlated noise which may reflect insufficient background potential drift removal.
Clearly, it is possible to retrieve similar and quantitative IP inversion models from integral chargeability data independent of current waveform by modelling the receiver and transmitter transfer functions. This means that the differences in the inversion models retrieved through standard inversions (Fig. 3) are mainly due to the lack of consideration of the current waveform, instead of to the fact Figure 6. Resistivity (top) and phase angle iCPA (second) inversion models with STDFs of 2 and 5 (white lines), misfit pseudo-sections for apparent resistivity (third) and IP (fourth) and apparent integral chargeability datum and forward response (bottom) for two responses marked in the pseudo-sections. Note the increase in integral chargeability misfit compared to the CPA and MPA parametrization (Figs 4 and 5) and the similarity between the phase angle models show in Fig. 4. Plot ranges and colour scales are the same as those used in Figs 4 and 5. Table 1. Summary of inversion model IP misfits for iCPA and corresponding misfits for CPA and MPA with gate STD propagated to integrated STD as uncorrelated (eq. 6) and correlated (eq. 7) error.  that the different pulse length excites different parts of the IP spectra. Furthermore, having phase angle as an IP model parameter can be beneficial since it enables direct comparison with data from lab IP measurements. These are carried out in the frequency domain and normally tabulate the phase angle (or the imaginary conductivity values). However, such inversions are limited to low-level parametrizations, which may be insufficient for explaining the data. In general, one way to mitigate the effects of these limitations is to consider the full-response IP in the inversion and invert for full response IP with whichever parametrization is suitable for explaining the data, for instance for CPA or MPA parametrizations.

C O N C L U S I O N
Time-domain induced polarization data is waveform dependant and data sets acquired with different current waveforms give different apparent chargeability data both for their integrals and for the full IP response. As demonstrated using field data with different pulse duration, the differences in the data propagate trough the inversion into the final models if the data is inverted without considering the current and potential waveforms. Thus, different IP inversion models are retrieved for different current waveforms and the retrieved IP models do not represent material properties of the subsurface. These can however still detect qualitative differences in IP magnitude.
To retrieve subsurface IP material properties quantitatively, we developed the iCPA inversion scheme, which models apparent integral chargeability data based on the CPA parametrization of the complex resistivity, expressed in terms of the medium resistivity and phase. The novel iCPA routine, apart from modelling the current and potential waveforms, also models the IP forward response integration. This modelling greatly reduces the bias in the inversion models due the current waveform. However, the iCPA inversion does not ensure that the CPA model is appropriate for describing the IP spectral content, and bias in the inversion models due to a non-appropriate parametrization of the IP phenomenon might still be present. The full-response inversion can be used for evaluating thoroughly the appropriateness of the iCPA inversion, also through the comparison of inversion schemes based on different spectral contents, such as the CPA and MPA inversions. Even so, the novel inversion scheme is a valuable method for mitigating effects of transmitter waveform when working with existing (older) data sets which lack full-response IP data and recognizing that full IP response data at present does not qualify as industry standard.

A C K N O W L E D G E M E N T S
The authors would like to thank Azadeh Rezvani for their enthusiastic support of this study, and especially the field survey.