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

Abstract

Summary We present a new algorithm for 3D forward modelling and spectral inversion of resistivity and time-domain full-decay induced polarization (IP) data. To our knowledge, all algorithms available for handling 3D 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}}$-${10^4}$ Hz) and then transform the response into the time domain, where we account for the transmitted current waveform. Solving Poisson's equation in 3D 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 (down-sampling) 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 down-sampling 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 3D anomalies and a field example, where lithology logs verify the results. The datasets 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 3D spectral inversion of resistivity and full-decay IP field data for both surface and borehole applications.

Integration of Electrical Resistivity Tomography and Induced Polarization for Characterization and Mapping of (Pb-Zn-Ag) Sulfide Deposits

“…This integration can enhance the interpretation of geophysical data and aid in various applications like groundwater exploration, mineral prospecting, and environmental studies. For full-decay IP modelling, the decay is divided into time gates and an apparent chargeability value is defined in each gate [22].…”

“…(a) the red dotted line represents the relationship between the applied current and time, while the blue line displays the behavior of the measured potential. (b)For full-decay IP modelling, the decay is divided into time gates and an apparent chargeability value is defined in each gate[22].…”

Three-dimensional time-lapse inversion of transient electromagnetic data, with application at an Icelandic geothermal site

“…The model nodes in the regular inversion mesh are defined on the mesh nodes with uniform node spacing in the horizontal direction and log-increasing node spacing in the vertical direction (Christensen et al 2017). The resistivity values in the cells of the forward meshes are obtained from the values of the inversion mesh through interpolation with an inverse distance weighting method (Madsen et al 2020). The log-increasing node spacing in the vertical direction helps balance the sensitivity of deep and shallow areas of the model to avoid favouring shallow TL changes, as would happen with linear vertical node spacing.…”

“…represents the Jacobian of the two acquisitions, translated from the forward mesh to the inversion mesh following Madsen et al (2020) and Zhang et al (2021); I is the identity matrix and R m and R u are the roughness matrices on model and model update. In eq.…”

PFLOTRAN-SIP: A PFLOTRAN Module for Simulating Spectral-Induced Polarization of Electrical Impedance Data

“…While there are numerous software for modeling geoelectrical data (e.g., IP-decay model [30], Res2Dinv [31][32][33], Aarhusinv [34], BERT [35,36], EarthImager3D [37], E4D [38], pyGIMLi [39], and ZondRes3D [40]), none capture the physics that are associated with dynamic subsurface processes. These software packages can also image frequency-dependent electrical conductivities, but they cannot capture dynamic subsurface processes.…”