ALBERT

Description: Prediction of stable mineral equilibria in the Earth's lithosphere is critical to unravel the tectonomagmatic history of exposed geological sections. While the recent advances in geodynamic modeling allow us to explore the dynamics of magmatic transfer in solid mediums, there is to date no available thermodynamic package that can easily be linked and efficiently be accounted for the computation of phase equilibrium in magmatic systems. Moreover, none of the existing tools fully exploit single point calculation parallelization, which strongly hinders their applicability for direct geodynamic coupling or for thermodynamic database inversions. Here, we present a new Mineral Assemblage Gibbs Energy Minimizer (magemin). The package is written as a parallel C library, provides a direct Julia interface, and is callable from any petrological/geodynamic tool. For a given set of pressure, temperature, and bulk‐rock composition magemin uses a combination of linear programming, extended Partitioning Gibbs Energy and gradient‐based local minimization to compute the stable mineral assemblage. We apply our new minimization package to the igneous thermodynamic data set of Holland et al. (2018), https://doi.org/10.1093/petrology/egy048 and produce several phase diagrams at supra‐solidus conditions. The phase diagrams are then directly benchmarked against thermocalc and exhibit very good agreement. The high scalability of magemin on parallel computing facilities opens new horizons, for example, for modeling reactive magma flow, for thermodynamic data set inversion, and for petrological/geophysical applications.

Description: Plain Language Summary: Understanding magmatic systems requires knowing how rocks melt. Because a single melting experiment can easily take weeks, it is impossible to do enough experiments to cover the whole range of pressure, temperature, and composition relevant for magmatic systems. We therefore need a way to interpolate in between conditions that are not directly covered by the experiments. It is long known that the best way to perform such interpolation is by using basic thermodynamic principles. For magmatic systems, this requires a well‐calibrated thermodynamic melting model. It also requires an efficient computational tool to predict the most stable configuration of minerals and melt. Since the 1980s, a number of such computational tools have been developed to perform a so‐called Gibbs energy minimization. These tools work very well for simpler systems but become very slow for recently developed, more realistic, melting models. Here, we describe a new method that combines some ideas of the previous methods with a new algorithm. Our method is faster and takes advantage of modern computer architectures. It can predict rock properties such as densities, seismic velocities, melt content, and chemistry. It can therefore be used to link physical observations with hard rock data of magmatic systems.

Description: Key Points: A new, parallel, Gibbs energy minimization approach is presented to compute multiphase multicomponent equilibria. It predicts parameters like stable phases, melt content, or seismic velocities as a function of chemistry and temperature/pressure conditions. Examples and benchmark cases are presented that apply the approach to magmatic systems.

Description: EC | H2020 | H2020 Priority Excellent Science | H2020 European Research Council (ERC) http://dx.doi.org/10.13039/100010663

Description: https://doi.org/10.5281/zenodo.6347567

Description: https://github.com/ComputationalThermodynamics/magemin.git

Description: Studies of host rock deformation around magmatic intrusions usually focus on the development of stresses directly related to the intrusion process. This is done either by considering an inflating region that represents the intruding body, or by considering multiphase deformation. Thermal processes, especially volume changes caused by thermal expansion are typically ignored. We show that thermal stresses around upper crustal magma bodies are likely to be significant and sufficient to create an extensive fracture network around the magma body by brittle yielding. At the same time, cooling induces decompression within the intrusion, which can promote the appearance of a volatile phase. Volatile phases and the development of a fracture network around the inclusion may thus be the processes that control magmatic‐hydrothermal alteration around intrusions. This suggests that thermal stresses likely play an important role in the development of magmatic systems. To quantify the magnitude of thermal stresses around cooling intrusions, we present a fully compressible 2D visco‐elasto‐plastic thermo‐mechanical numerical model. We utilize a finite difference staggered grid discretization and a graphics processing unit based pseudo‐transient solver. First, we present purely thermo‐elastic solutions, then we include the effects of viscous relaxation and plastic yielding. The dominant deformation mechanism in our models is determined in a self‐consistent manner, by taking into account stress, pressure, and temperature conditions. Using experimentally determined flow laws, the resulting thermal stresses can be comparable to or even exceed the confining pressure. This suggests that thermal stresses alone could result in the development of a fracture network around magmatic bodies.

Description: Plain Language Summary: Quantifying the stresses that magma bodies exert on the surrounding rocks is an important part of understanding mechanical processes that control the evolution of magmatic systems and volcanic eruptions. Previous analytical or numerical models typically describe the mechanical response to changes in magma volume due to intrusion or extraction of magma. However, volume changes related to thermal expansion/contraction around a cooling magma body are often neglected. Here, we develop a new software which runs on modern graphics processing unit machines, to quantity the effect of this process. The results show that stresses due to thermal expansion/contraction are significant, and often large enough to fracture the rocks nearby the magma body. Such fracture networks may form permeable pathways for the magma or for fluids such as water and CO〈sub〉2〈/sub〉, thus influencing the evolution of magmatic and hydrothermal systems. Finally we show that cooling and shrinking of magma bodies causes significant decompression which can influence the chemical evolution of the magma during crystallization and devolatilization.

Description: Key Points: We present a numerical quantification of the effect of thermal stresses in visco‐elasto‐plastic rock with tensile and dilatant shear failure. The pressure drop in thermally contracting upper crustal magma bodies can exceed 100 MPa, potentially triggering devolatilization. Thermal cracking can create an extensive fracture network around an upper crustal magma body.

Description: European Research Council http://dx.doi.org/10.13039/501100000781

Description: https://zenodo.org/record/6958273

Description: https://doi.org/10.5281/zenodo.6958273

Description: One of the main methods to determine the strength of the lithosphere is by estimating it's effective elastic thickness. This method assumes that the lithosphere is a thin elastic plate that floats on the mantle and uses both topography and gravity anomalies to estimate the plate thickness. Whereas this seems to work well for oceanic plates, it has given controversial results in continental collision zones. For most of these locations, additional geophysical data sets such as receiver functions and seismic tomography exist that constrain the geometry of the lithosphere and often show that it is rather complex. Yet, lithospheric geometry by itself is insufficient to understand the dynamics of the lithosphere as this also requires knowledge of the rheology of the lithosphere. Laboratory experiments suggest that rocks deform in a viscous manner if temperatures are high and stresses low, or in a plastic/brittle manner if the yield stress is exceeded. Yet, the experimental results show significant variability between various rock types and there are large uncertainties in extrapolating laboratory values to nature, which leaves room for speculation. An independent method is thus required to better understand the rheology and dynamics of the lithosphere in collision zones. The goal of this paper is to discuss such an approach. Our method relies on performing numerical thermomechanical forward models of the present-day lithosphere with an initial geometry that is constructed from geophysical data sets. We employ experimentally determined creep-laws for the various parts of the lithosphere, but assume that the parameters of these creep-laws as well as the temperature structure of the lithosphere are uncertain. This is used as a priori information to formulate a Bayesian inverse problem that employs topography, gravity, horizontal and vertical surface velocities to invert for the unknown material parameters and temperature structure. In order to test the general methodology, we first perform a geodynamic inversion of a synthetic forward model of intraoceanic subduction with known parameters. This requires solving an inverse problem with 14–16 parameters, depending on whether temperature is assumed to be known or not. With the help of a massively parallel direct-search combined with a Markov Chain Monte Carlo method, solving the inverse problem becomes feasible. Results show that the rheological parameters and particularly the effective viscosity structure of the lithosphere can be reconstructed in a probabilistic sense. This also holds, with somewhat larger uncertainties, for the case where the temperature distribution is parametrized. Finally, we apply the method to a cross-section of the India–Asia collision system. In this case, the number of parameters is larger, which requires solving around 1.9  x 10 6 forward models. The resulting models fit the data within their respective uncertainty bounds, and show that the Indian mantle lithosphere must have a high viscosity. Results for the Tibetan plateau are less clear, and both models with a weak Asian mantle lithosphere and with a weak Asian lower crust fit the data nearly equally well.

Description: The mechanisms that result in the formation of high-pressure (HP) and ultrahigh-pressure (UHP) rocks are controversial. The usual interpretation assumes that pressure is close to lithostatic, petrological pressure estimates can be transferred to depth, and (U)HP rocks have been exhumed from great depth. An alternative explanation is that pressure can be larger than lithostatic, particularly in continental collision zones, and (U)HP rocks could thus have formed at shallower depths. To better understand the mechanical feasibility of these hypotheses, we performed thermomechanical numerical simulations of a typical subduction and collision scenario. If the subducting crust is laterally homogeneous and has small effective friction angles (and is thus weak), we reproduce earlier findings that 〈20% deviation of lithostatic pressure occurs within a subduction channel. However, many orogenies involve rocks that are dry and strong, and the crust is mechanically heterogeneous. If these factors are taken into account, simulations show that pressures can be significantly larger than lithostatic within nappe-size, mechanically strong crustal units, or within a strong lower crust, as a result of tectonic deformation. Systematic simulations show that these effects are most pronounced at the base of the crust (at ~40 km), where pressures can reach 2–3 GPa (therefore within the coesite stability field) for millions of years. These pressures are often released rapidly during ongoing deformation. Relating metamorphic pressure estimates to depth might thus be problematic in mechanically heterogeneous crustal rock units that appear to have been exhumed in an ultrafast manner.

Description: Many unresolved questions in geodynamics revolve around the physical behaviour of the two-phase system of a silicate melt percolating through and interacting with a tectonically deforming host rock. Well-accepted equations exist to describe the physics of such systems and several previous studies have successfully implemented various forms of these equations in numerical models. To date, most such models of magma dynamics have focused on mantle flow problems and therefore employed viscous creep rheologies suitable to describe the deformation properties of mantle rock under high temperatures and pressures. However, the use of such rheologies is not appropriate to model melt extraction above the lithosphere–asthenosphere boundary, where the mode of deformation of the host rock transitions from ductile viscous to brittle elasto-plastic. Here, we introduce a novel approach to numerically model magma dynamics, focusing on the conceptual study of melt extraction from an asthenospheric source of partial melt through the overlying lithosphere and crust. To this end, we introduce an adapted set of two-phase flow equations, coupled to a visco-elasto-plastic rheology for both shear and compaction deformation of the host rock in interaction with the melt phase. We describe in detail how to implement this physical model into a finite-element code, and then proceed to evaluate the functionality and potential of this methodology using a series of conceptual model setups, which demonstrate the modes of melt extraction occurring around the rheological transition from ductile to brittle host rocks. The models suggest that three principal regimes of melt extraction emerge: viscous diapirism, viscoplastic decompaction channels and elasto-plastic dyking. Thus, our model of magma dynamics interacting with active tectonics of the lithosphere and crust provides a novel framework to further investigate magmato-tectonic processes such as the formation and geometry of magma chambers and conduits, as well as the emplacement of plutonic rock complexes.

Description: Numerical models of mantle convection typically employ a temperature- or pressure-dependent viscous or viscoplastic rheology and a free slip upper boundary condition. The Earth, however, has a stress-free rather than a free slip surface condition. In addition, with decreasing temperature, the viscosity of rocks increases, which might induce a change from viscous to elastic behaviour (depending on the timescale of deformation). Here, we study the effects of both a Maxwell viscoelastic rheology and a free surface upper boundary condition on viscoelastic convection with a strongly temperature dependent rheology. We particularly focus on the effect of elasticity on the stress state of the lithosphere. Results show that convection vigor and heat transport are not significantly altered by the upper boundary condition or by elasticity. However, the stress state of the lithosphere is significantly affected by both factors. If elasticity is unimportant, a free surface upper boundary condition results in significantly elevated surface stresses (which are up to two magnitudes larger than in the free slip case). Elasticity counteracts this effect and significantly reduces the surface stresses, but distributes stresses over a thicker layer than in the case of a purely viscous rheology. At Earth-like conditions, this effect is significant. While it is warranted to use a free slip upper boundary condition and neglecting elasticity when studying mantle convection and its effect on the thermal state of the Earth, both factors are significant when one wants to predict the stress state of the lithosphere and related questions. Additional 2-D simulations of a plume impinging on a constant thickness and constant viscosity lithosphere show that reasonable parameters might induce lithospheric stress levels that are on the order of a GPa or larger, for viscous free surface models, and that these stresses are several orders of magnitude larger than stresses that occur for free slip models. This suggests that the effect of a free surface should not be ignored in models where the rheology is stress-dependent, such as in viscoplastic models of self-consistent plate tectonics.

Description: Many salt diapirs are thought to have formed as a result of down-building, which implies that the top of the diapir remained close to the surface during syn-halokinetic sediment deposition. Down-building is largely a 3-D process and in order to better understand what controls the patterns of the diapirs that form by this process, we here perform 3-D numerical models of down-built diapirs initiated by the gravity instability in linear viscous materials and compare the results with analytical models. We vary several parameters of the numerical models such as initial salt thickness, sedimentation rate, salt viscosity, salt-sediment viscosity ratio as well as the density of sediments. Down-building of 3-D diapirs only occurs for a certain range of parameters and is favoured by lower sediment/salt viscosity contrasts and sedimentation rates in agreement with analytical predictions and findings from previous 2-D models. However, the models show that the sedimentation rate has an additional effect on the formation and evolution of 3-D diapir patterns. At low sedimentation rates, salt ridges that form during early model stages remain preserved at later stages as well. For higher sedimentation rates, the initial salt ridges are covered up and finger-like diapirs form at their junctions, which results in different salt exposure patterns at the surface. Once the initial pattern of diapirs is formed, higher sedimentation rate can also result in covered diapirs if the diapir extrusion velocity is insufficiently large. We quantify the effect of sedimentation rate on the number of diapirs exposed at the surface as well as on their spacing and we explain the observations with analytical predictions using thick-plate analytical models. In some cases, this final pattern is distinctly different from the initial polygonal pattern.

Description: SUMMARY Imaging subsurface structures, such as salt domes, magma reservoirs or subducting plates, is a major challenge in geophysics. Seismic imaging methods are, so far, the most precise methods to open a window into the Earth. However, the methods may not yield the exact depth or size of the imaged feature and may become distorted by phenomena such as seismic anisotropy, fluid flow, or compositional variations. A useful complementary method is therefore to simulate the mechanical behaviour of rocks on large timescales, and compare model predictions with observations. Recent studies have used the (non-linear) Stokes equations and geometries from seismic studies in combination with an adjoint-based approach to invert for rheological parameters that are consistent with surface observations such as GPS velocities. Nevertheless, it would be useful to use other surface observations, such as principal stress directions, as constraints as well. Here, we derive the adjoint formulation for the case that principal stress directions are used as observables with respect to rheological parameters. Both an algebraic and a discretized derivation of the adjoint equations are described. This thus enables the usage of two data fields - surface velocities and stress directions - as a misfit for the inversion. We test the performance of the inversion for principal stress directions on simplified 3-D test cases. Finally, we demonstrate how the adjoint approach can be used to compute 3-D geodynamic sensitivity kernels, which highlight the areas in the model domain that have the largest impact on the misfit value of a particular point. This provides a simple, yet powerful, way to visualize which parts of the model domain are of key importance if changing rheological constants.