2D Hamiltonian Monte Carlo joint inverse modeling of potential field data

Here we present preliminary results from an algorithm for joint inverse modeling of gravity and total magnetic anomaly (TMA) data aiming to target the density distribution and the magnetization in the subsurface of the Earth by means of the Hamiltonian Monte Carlo method (HMC) (Neal, 2011; Fichtner et al., 2018). The model parameterization is defined in terms of 2D polygonal bodies characterized by uniform density and magnetizations (induced and/or remnant) and the gravity and magnetic non-linear forward calculations are performed using the formulas of Talwani et al. (1959) and Talwani & Heirtzler (1962) respectively, the latter ones checked in detail recently by Ghirotto et al. (2021). The main benefits using polygonal bodies are i) a huge reduction of the model space size for well-defined geological bodies compared to gridded approaches and ii) a computationally much cheaper algorithm in terms of both memory requirements and number of calculations. The unknown parameters are represented by both the positions of the polygon vertices and their density and magnetizations, with the option of limiting the inversion to any of them. Following the probabilistic approach to inverse problems, the goal of the HMC inversion strategy is to explore the posterior probability density of the model parameters (PPD), obtaining as a result a collection of models representing samples of the PPD. In addition, statistical analysis performed on this collection could provide useful measures of parameters uncertainty and plausible geophysical scenarios, which cannot be handled using traditional deterministic inversion methods. To help steering the inversion process toward high-probability areas in the model space manifold, HMC requires the computation of the gradient of the PPD with respect to the model parameters, achieved here by employing the technique of automatic differentiation.