Hamiltonian Monte Carlo applied to inverse petrological problems

The petrogenesis of rocks can be investigated by integrating multiple datasets and numerical methods. Petrology benefits from the abundance of field, petrographic, geochronological and geochemical data. Moreover, the growing number of numerical methods opens new opportunities for inversion. Inversion allows the quantification of useful parameters, offering deeper insights into natural processes. However, an increasing number of parameters to invert for is accompanied by significant computational costs. Therefore, relevant algorithms are needed to perform inversion and uncertainty quantification in high-dimensional spaces. Here, we demonstrate the use of Hamiltonian Monte Carlo for inverse diffusion modeling in a petrological framework. Hamiltonian Monte Carlo is a gradient-based method that efficiently explores high-dimensional parameter spaces. We implemented it using the Turing.jl package in Julia which makes use of Automatic Differentiation to efficiently explore the parameter space. Our analysis focused on the calculation of the initial cooling rate, equilibration temperature and effective mineral grain size. We fit compositional garnet data and Ar-muscovite geochronological data from the Pindos metamorphic sole (Greece) by using two different forward diffusion models. Our joint inversion shows an initial equilibration temperature of 632.3 ± 9.3 °C and a cooling rate of 202.6 ± 72.0 °C/Myr. These values reproduce not only thermobarometric and geochronological observations but also fit the garnet composition profiles and the ⁴⁰Ar/³⁹Ar age of muscovite. We finally aim to highlight the potential of Hamiltonian Monte Carlo as a robust method to perform high-dimensional inversion and constrain complex petrological processes.