Potential fields data modeling: new frontiers in forward and inverse problems

Since the '50, potential fields data modeling has played an important role in analyzing the density and magnetization distribution in Earth's subsurface for a wide variety of applications. Examples are the characterization of ore deposits and the assessment of geothermal and petroleum potential, which turned out to be key contributors for the economic and industrial development after World War II. Current modeling methods mainly rely on two popular parameterization approaches, either involving a discretization of target geological bodies by means of 2D to 2.75D horizontal prisms with polygonal vertical cross-section (polygon-based approach) or prismatic cells (prism-based approach). Despite the great endeavour made by scientists in recent decades, inversion methods based on these parameterization approaches still suffers from a limited ability to (i) realistically characterize the variability of density and magnetization expected in a study area and (ii) take into account the strong non-uniqueness affecting potential fields theory. The prism-based approach is used in linear deterministic inverse methods, which provide just one single solution, preventing uncertainty estimation and statistical analysis on the parameters we would like to characterize (i.e, density or magnetization). On the contrary, the polygon-based approach is almost exclusively exploited in a trial-and-error modeling strategy, leaving the potential to develop innovative inverse methods untapped. The reason is two-fold, namely (i) its strongly non-linear forward problem requires an efficient probabilistic inverse modeling methodology to solve the related inverse problem, and (ii) unpredictable cross-intersections between polygonal bodies during inversion represent a challenging task to be tackled in order to achieve geologically plausible model solutions. The goal of this thesis is then to contribute to solving the critical issues outlined above, developing probabilistic inversion methodologies based on the polygon- and prism-based parameterization approaches aiming to help improving our capability to unravel the structure of the subsurface. Regarding the polygon-based parameterization strategy, at first a deep review of its mathematical framework has been performed, allowing us (i) to restore the validity of a recently criticized mathematical formulation for the 2D magnetic case, and (ii) to find an error sign in the derivation for the 2.75D magnetic case causing potentially wrong numerical results. Such preliminary phase allowed us to develop a methodology to independently or jointly invert gravity and magnetic data exploiting the Hamilton Monte Carlo approach, thanks to which collection of models allow researchers to appraise different geological scenarios and fully characterize uncertainties on the model parameters. Geological plausibility of results is ensured by automatic checks on the geometries of modelled bodies, which avoid unrealistic cross-intersections among them. Regarding the prism-based parameterization approach, the linear inversion method based on the probabilistic approach considers a discretization of target geological scenarios by prismatic bodies, arranged horizontally to cover it and finitely extended in the vertical direction, particularly suitable to model density and magnetization variability inside strata. Its strengths have been proven, for the magnetic case, in the characterization of the magnetization variability expected for the shallower volcanic unit of the Mt. Melbourne Volcanic Field (Northern Victoria Land, Antarctica), helping significantly us to unravel its poorly known inner geophysical architecture

HMCLab: a framework for solving diverse geophysical inverse problems using the Hamiltonian Monte Carlo method

The use of the probabilistic approach to solve inverse problems is becoming more popular in the geophysical community, thanks to its ability to address nonlinear forward problems and to provide uncertainty quantification. However, such strategy is often tailored to specific applications and therefore there is a lack of a common platform for solving a range of different geophysical inverse problems and showing potential and pitfalls. We demonstrate a common framework to solve such inverse problems ranging from, e.g, earthquake source location to potential field data inversion and seismic tomography. Within this approach, we can provide probabilities related to certain properties or structures of the subsurface. Thanks to its ability to address high-dimensional problems, the Hamiltonian Monte Carlo (HMC) algorithm has emerged as the state-of-the-art tool for solving geophysical inverse problems within the probabilistic framework. HMC requires the computation of gradients, which can be obtained by adjoint methods, making the solution of tomographic problems ultimately feasible. These results can be obtained with "HMCLab", a tool for solving a range of different geophysical inverse problems using sampling methods, focusing in particular on the HMC algorithm. HMCLab consists of a set of samplers and a set of geophysical forward problems. For each problem its misfit function and gradient computation are provided and, in addition, a set of prior models can be combined to inject additional information into the inverse problem. This allows users to experiment with probabilistic inverse problems and also address real-world studies. We show how to solve a selected set of problems within this framework using variants of the HMC algorithm and analyze the results. HMCLab is provided as an open source package written both in Python and Julia, welcoming contributions from the community.Comment: 21 pages, 4 figure

Geophysical imaging unveils the largest pull-apart basin in East Antarctica

West Antarctica hosts one of the largest continental rift systems on Earth, the West Antarctic Rift System (WARS) that forms the lithospheric cradle for the West Antarctic Ice Sheet. The WARS is known to have experienced several stages of extension starting with distributed/wide mode extension in the Cretaceous, followed by narrower mode and variably oblique extension in the Cenozoic, the latter potentially triggered by the onset of oceanic seafloor spreading in the Adare Basin (Davey et al., 2016, GRL). However, the extent and impact of Cenozoic extension and transtension within the Transantarctic Mountains sector of East Antarctica is much less well understood. Here we present results from a new project (REGGAE) that by analysing aeromagnetic, aerogravity and land-gravity and bedrock topography images and models provides key new geophysical constraints on the form, extent and kinematics of the largest Cenozoic pull-apart basin recognised so far in East Antarctica, the Rennick Graben (RG). Potential field imaging reveals the extent of part of a Jurassic tholeiitic Large Igneous Province preserved within the RG and helps delineate the inherited structural architecture of the underlying Ross-age basement in northern Victoria Land, including highly magnetic arc basement in the northern Wilson Terrane and the subglacial extent of a thrust fault belt located between the western flank of the RG and the eastern margin of Wilkes Subglacial Basin (WSB). We show that the RG is a major composite right-lateral pull-part basin that extends from the Oates Coast to the Southern Cross Mountains crustal block and propose that it is kinematically connected with both the western edge of the WARS and the eastern margin of the WSB. More cryptic evidence for an earlier phase of left-lateral strike slip deformation is also emerging from our recent geological field work in the study region and relatively subtle offsets in aeromagnetic anomaly patterns. Our findings suggest that the RG is part of a distributed region of the continental lithosphere in East Antarctica that was preferentially deformed in response to Cenozoic transtensional stresses that likely also facilitated propagation of accelerated oceanic transform faulting in the adjacent oceanic lithosphere located between southeastern Australia and Tasmania

Enhanced images and new models of the Wilkes Subglacial Basin help constrain the variability in geological boundary conditions for the East Antarctic Ice Sheet

The Wilkes Subglacial Basin (WSB) is a huge tectonic feature formed by Cenozoic lithospheric flexure coupled with Mesozoic to Cenozoic extension localised in sub-basins (Paxman et al., 2019, JGR). The deep northern WSB underlies the catchments of the Matusevich, Cook, Ninnis and Mertz glaciers that are largely marine-based, which renders them more vulnerable to past and predicted future ocean and climate warming. Here we present airborne radar and enhanced magnetic and gravity views of the northern WSB that help unveil the spatial variability in geological boundary conditions for this key sector of the East Antarctic Ice Sheet (EAIS). Residual gravity anomalies obtained by stripping out Moho effects were compared with aeromagnetic anomaly images to glean new perspectives into intra-crustal features. Depth to magnetic and gravity source estimates were then used to help derive the first combined 2D forward models for the region. We first examine a model crossing the northern WSB extending from the Matusevich Glacier to the deep Cook Basins. The model reveals a major crustal boundary along the eastern margin of the WSB interpreted as separating the Ross Orogen from a composite Precambrian Wilkes Terrane buried beneath Devonian to Jurassic sediments and early Cambrian metasediments. By analogy with the better understood Rennick Graben in northern Victoria Land, the Cook basins are interpreted as glacially over deepened grabens. The Cook basins clearly play a major role in EAIS dynamics, as they steer fast glacial flow deep into the interior of East Antarctica where they connect to the Central Basins. Our new model across these basins shows that the inferred Precambrian basement is both shallower and of more felsic bulk composition compared to the Cook basins. This fundamental difference in basement depth, bulk composition and thickness of sedimentary cover is likely to exert major influences on geothermal heat variability in this key sector of the EAIS

Focused Inversion of Gravimetric and Magnetotelluric Data for Geothermal Investigations

Focused inversion techniques may be applied to geophysical data inversion in order to image complex structures in the subsoil. These algorithms may image complex \u201cblocky\u201d structures giving useful information in geothermal exploration that may be smoothed out by standard inversion algorithms that use stabilizer that penalize sharp transitions. We have tested the modified total variation and the maximum gradient support stabilizers in the inversion of synthetic and field magnetotelluric and gravimetric data. The gravimetric data from the Luhoi geothermal prospect have been used to map the sharp density transition between the sandstone and the overlying claystone layers. The resulting horst structure imaged in 2D and 3D models by the maximum gradient support stabilizer solution allow to trace the main fault system that drives the up-flow of hydrothermal waters. The 1D magnetotelluric \u201cblocky\u201d models with lateral constrain (pseudo-3D) image the lithological contact between the claystone and sandstone far from the horst area and reveal resistivity variations in the claystone layer associated with sand lenses. In the horst area, resistivity models image hydrothermal alteration affecting the sandstone layer

HMCLab: a framework for solving diverse geophysical inverse problems using the Hamiltonian Monte Carlo method

The use of the probabilistic approach to solve inverse problems is becoming more popular in the geophysical community, thanks to its ability to address nonlinear forward problems and to provide uncertainty quantification. However, such strategy is often tailored to specific applications and therefore there is a need for common platforms to solve different geophysical inverse problems and showing potential and pitfalls of the methodology. In this work, we demonstrate a common framework within which it is possible to solve such inverse problems ranging from, for example, earthquake source location to potential field data inversion and seismic tomography. This allows us to fully address nonlinear problems and to derive useful information about the subsurface, including uncertainty estimation. This approach can, in fact, provide probabilities related to certain properties or structures of the subsurface, such as histograms of the value of some physical property, the expected volume of buried geological bodies or the probability of having boundaries defining different layers. Thanks to its ability to address high-dimensional problems, the Hamiltonian Monte Carlo (HMC) algorithm has emerged as the state-of-the-art tool for solving geophysical inverse problems within the probabilistic framework. HMC requires the computation of gradients, which can be obtained by adjoint methods. This unique combination of HMC and adjoint methods is what makes the solution of tomographic problems ultimately feasible. These results can be obtained with ‘HMCLab’, a numerical laboratory for solving a range of different geophysical inverse problems using sampling methods, focusing in particular on the HMC algorithm. HMCLab consists of a set of samplers (HMC and others) and a set of geophysical forward problems. For each problem its misfit function and gradient computation are provided and, in addition, a set of prior models can be combined to inject additional information into the inverse problem. This allows users to experiment with probabilistic inverse problems and also address real-world studies. We show how to solve a selected set of problems within this framework using variants of the HMC algorithm and analyse the results. HMCLab is provided as an open source package written both in Python and Julia, welcoming contributions from the community.ISSN:0956-540XISSN:1365-246

Hamiltonian Monte Carlo Probabilistic Joint Inversion of 2D (2.75D) Gravity and Magnetic Data

Two-dimensional modeling of gravity and magnetic anomalies in terms of polygonal bodies is a popular approach to infer possible configurations of geological structures in the subsurface. Alternatively to the traditional trial-and-error manual fit of measured data, here we illustrate a probabilistic strategy to solve the inverse problem. First we derive a set of formulae for solving a 2.75-dimensional forward model, where the polygonal bodies have a given finite lateral extent, and then we devise a Hamiltonian Monte Carlo algorithm to jointly invert gravity and magnetic data for the geometry and properties of the polygonal bodies. This probabilistic approach fully addresses the nonlinearity of the forward model and provides uncertainty estimation. The result of the inversion is a collection of models which represent the posterior distribution, analysis of which provides estimates of sought properties and may reveal different scenarios.ISSN:0094-8276ISSN:1944-800

Magnetic Anomalies Caused by 2D Polygonal Structures With Uniform Arbitrary Polarization: New Insights From Analytical/Numerical Comparison Among Available Algorithm Formulations

Since the '60s of the last century, the calculation of the magnetic anomalies caused by 2D uniformly polarized bodies with polygonal cross\u2010section has been mainly performed using the popular algorithm of Talwani and Heirtzler (1962, 1964). Recently, Kravchinsky et al. (2019, https://doi.org/10.1029/2019GL082767) claimed errors in the above algorithm formulation, proposing new corrective formulas and questioning the effectiveness of almost 60 years of magnetic calculations. Here we show that the two approaches are equivalent and Kravchinsky et al.'s formulas simply represent an algebraic variant of those of Talwani and Heirtzler. Moreover, we analyze a large amount of random magnetic scenarios, involving both changing\u2010shape polygons and a realistic geological model, showing a complete agreement among the magnetic responses of the two discussed algorithms and the one proposed by Won and Bevis (1987, https://doi.org/10.1190/1.1442298). We release the source code of the algorithms in Julia and Python languages

Imaging the Mt. Melbourne Volcanic Field (Northern Victoria Land, Antarctica): a Hamiltonian Monte Carlo approach applied to high-resolution aeromagnetic data

The Mt. Melbourne Volcanic Complex (MMVC) is located in Northern Victoria Land (Antarctica) along the western flank of the West Antarctic Rift System, at the boundary with the Transantarctic Mountains. It is constituted by two main volcanic areas, i.e. the Mt. Melbourne Edifice (MME) and the Cape Washington Shield (CWS), and some other minor centres. To date, the inner structure of this volcanic complex is still poorly known, being the direct geological information on site confined to either glacial erratics or a few rock outcrops not hidden by the ice sheet. Consequently, even the temporal building up and evolution of the MMVC as well as its primary magmatic source are still under investigation (debated). Recently, we attempted to define the geological structure of the MMVC by means of digital enhancement and forward modeling performed on a high-resolution aeromagnetic dataset (Ghirotto et al. 2020, EGU). Coupling both information derived from past geological/geophysical studies and unpublished magnetic susceptibility measurements from rock samples collected in the field, we proposed two models to explain the chronological evolution of the MME and CWS. These models involve either i) major magmatic events occurred in periods of both normal and reverse magnetic polarity or ii) only magmatic flows with normal polarity. To gain further insights into the geological structure and the geodynamic evolution of the MMVC in relation to the two proposed models, we develop here a Hamiltonian Monte Carlo (HMC) algorithm (Fichtner et al. 2018) based on the probabilistic approach to inverse problems. To date, this methodology has never been applied to aeromagnetic data for geological studies. In detail, the above proposed models provide some soft a priori information from which to start exploring potential solutions. The parameterization of the volcanic area is defined in terms of 2-D polygonal bodies, representing e.g. magmatic lava flows, where the unknown parameters are represented by both the position of the vertices and/or the magnetization (induced and/or remnant), resulting in a non-linear forward model. The HMC algorithm requires the computation of gradients of the posterior probability density (PPD), i.e., derivatives of the objective functional with respect to the position of vertices of the bodies and magnetization, in order to better move the inversion process toward high-probability areas in the model space manifold. We implement such calculations using automatic differentiation, a tool which is very accurate and fast compared to other approaches such as finite difference. The result of the inversion is then a collection of models representing the PPD, from which statistical analysis can provide measures of uncertainty and plausible geological scenarios. In this study we present some preliminary results of applying the above-mentioned methodology, which finally could help unravel the framework of the MMVC