A new and efficient computational thermodynamics approach for magmatic systems

During the last decade, the development of numerical geodynamic tools helped the geosciences community to unravel complex thermo-mechanical processes at play during plate tectonics. Yet, the high computational cost of thermodynamic calculations, which simulates phase change, hampers our ability to integrate complex chemistry in such problems. This is particularly important for simulating magmatic processes, where the chemistry of differentiating melts can vary significantly from the mantle to the upper crust. The typical approach, currently used, is to precompute one or many phase diagrams and use them as look-up tables. For many geodynamic processes this is adequate but when the melt chemistry varies drastically it would be better to be able to do thermodynamic calculations on the fly, along with the geodynamic models.

For that, the thermodynamic computational approach must be sufficiently fast, should work fully automatically and be tuned for melting models of magmatic systems, for example by utilizing the recently developed thermodynamic melting model of Holland et al. (2018). Existing approaches do not fulfill all criteria, which is why we have developed a new computational library for this purpose. Our code is written in C, runs on massively parallel machines (MPI) and uses an adaptive mesh refinement strategy to compute phase diagrams. At the moment we have focused on the 'igneous set' of the Holland & Powell dataset (as defined in the thermocalc software) to calculate stable phase equilibria in the system K2O–Na2O–CaO–FeO–MgO–Al2O3–SiO2–H2O–TiO2–Fe2O3–Cr2O3 (KNCFMASHTOCr). The code uses pressure, temperature and bulk-rock composition as input and returns relevant petrological and geodynamic information such as (but not restricted to) stable assemblage, phase fractions and phase densities. Different than many of the existing approaches, our method can efficiently utilize initial guesses which naturally occur in geodynamic simulations where the changes in chemistry between timesteps are usually minor.

The methodology performs a Gibbs free energy minimization and involves two main steps. First, we use a combination of levelling methods (iterative change of base) to reduce the number of potential (pure and solution) phases and to bring the G-hyperplane close to solution. Second, we use a partitioning of Gibbs energy approach coupled with local minimization to satisfy the Gibbs-Duhem rule and to retrieve the final set of stable solution phases. To illustrate the efficiency of the library up to supra-solidus conditions we present a set of dry phase diagrams and compare results of our computations with thermocalc calculations.

Ongoing development includes the treatment of solvus to extend its applicability to complex wet systems involving solution phase such as amphibole.