Towards fully-coupled thermodynamic-thermomechanical two-phase flow models of transcrustal magmatic systems

The modern view of magmatic systems includes transport and storage of melt at depths within the solid crust. Parameters that control the depths and duration of magma storage include the rheology (especially the depth of the brittle-ductile transition) and chemistry of the melt and solid, the crustal heat budget, and the tectonic setting. Understanding the large spatio-temporal variation in the magmatic output of volcanic arcs, therefore, requires consideration of solid deformation, melt transport, and chemical differentiation in transcrustal magmatic systems. We present preliminary reactive transport models of transcrustal magmatic systems in the ductile limit. The model couples deformation of a compressible solid and Darcy-flow of melt through the solid matrix to phase equilibria calculations. The thermomechanical solver is self-developed in the Julia language and employs a matrix-free pseudo transient solution technique (e.g., Raess et al. 2022). Julia solves the two-language problem and allows for rapid transition from the development stage to massively parallelized production simulations. The thermomechanical two-phase flow equations are numerically discretized on 1D and 2D cartesian finite difference grids. Two different approaches of thermodynamic-thermomechanical coupling are implemented: (i) On-the-fly calculation of stable mineral phases using the Gibbs energy minimization software MAGEMin (Riel et al. 2022) and (ii) interpolation from a precompiled phase diagram. MAGEMin is parallelized and shipped with a Julia wrapper. These features allow near seamless integration of thermodynamic calculations into the newly-developed geodynamic algorithm. The following workflow is applied: (1) Stable mineral phases are calculated at each numerical grid point as a function of local pressure and temperature values using MAGEMin. A starting bulk rock composition characteristic to arc magmas is used for the first time step. (2) The densities as well as the predicted major-oxide fractions of the solid and melt are used in the two-phase flow algorithm to solve for solid deformation, porosity evolution, and major-oxide transport. (3) The local major-oxide fractions predicted by the thermomechanical solver are used as new starting compositions in MAGEMin for the next time step.  For the second approach, we compiled data from published fractional crystallization experiments at relevant pressure and temperature conditions. The data set contains the major-oxide fractions as well as calculated melt and solid densities. During the thermomechanical solution procedure, the data is interpolated from this phase diagram according to local pressure and temperature conditions at each grid point. A detailed comparison of the two approaches and potential implications for natural systems are presented.  References:  Räss, L., Utkin, I., Duretz, T., Omlin, S., & Podladchikov, Y. Y. (2022). Assessing the robustness and scalability of the accelerated pseudo-transient method. Geoscientific Model Development, 15(14), 5757-5786. Riel, N., Kaus, B. J., Green, E. C. R., & Berlie, N. (2022). MAGEMin, an efficient Gibbs energy minimizer: application to igneous systems. Geochemistry, Geophysics, Geosystems, 23(7), e2022GC010427.