3-D time-domain spectral inversion of resistivity and full-decay induced polarization data - full solution of Poisson's equation and modelling of the current waveform

We present a new algorithm for 3-D forward modelling and spectral inversion of resistivity and time-domain full-decay induced polarization (IP) data. To our knowledge, all algorithms available for handling 3-D spectral inversion of full-decay IP data use a time-domain approximation to Poisson's equation in the forward response. To avoid this approximation, we compute the response in the frequency domain solving the full version of Poisson's equation for a range of frequencies (10-8-104 Hz) and then transform the response into the time domain, where we account for the transmitted current waveform. Solving Poisson's equation in 3-D is computationally expensive and in order to balance accuracy, time, and memory usage we introduce the following: (1) We use two separate meshes for the forward response and the model update, respectively. The forward mesh is an unstructured tetrahedral mesh allowing for local refinements whereas the model (inversion) mesh is a node-based structured mesh, where roughness constraints are easily implemented. By decoupling the two meshes, they can be tuned for optimizing the forward accuracy and the inversion resolution, independently. (2) A singularity removal method known from resistivity modelling has been adapted to the complex IP case and is applied to minimize the numerical errors caused by the fast changing potential close to the source electrodes. The method includes splitting the potential field into a primary part (response of a homogenous background) and a secondary part (from the anomalies). Two different forward meshes are then used to compute the forward response: a dense mesh for the primary potential field (only computed once for each frequency) and a coarser mesh for the secondary potential field (computed in each iteration step of the inversion). With this method, the singularity is minimized and the memory usages is decreased significantly at the same time. (3) Finally, we are sparsing (downsampling) the Jacobian matrix based on a threshold value of the normalized sensitivity. The Jacobian computation is performed by time-transforming the frequency-domain Jacobian obtained through the adjoint method. The Jacobian downsampling is carried out before the time-transform in the frequency domain, thus avoiding the time-transformation of the Jacobian elements with negligible sensitivity. We invert resistivity data and all IP time-gates simultaneously and use the Gauss-Newton model update to minimize the L2 misfit function. We invert the resistivity data and all IP time-gates simultaneously and use the Gauss-Newton model update to minimize the L2 misfit function. We demonstrate the performance of our inversion approach with a synthetic data example with 3-D anomalies and a field example, where lithology logs verify the results. The data sets contain 1256 quadrupole measurements with 33 IP time-gates each. The inversions results show good data fits and model retrieval. The inversion takes approximately one hour per iteration using four CPUs. With this speed and accuracy, we believe this modelling and inversion approach will be a strong tool for 3-D spectral inversion of resistivity and full-decay IP field data for both surface and borehole applications.