Diapirism is crucial for heat and mass transfer in many geodynamic processes. Understanding diapir ascent velocity is vital for assessing its significance in various geodynamic settings. Although analytical estimates exist for ascent velocities of diapirs in power-law viscous, stress weakening fluids, they lack validation through 3D numerical calculations. Here, we improve these estimates by incorporating combined linear and power-law viscous flow and validate them using 3D numerical calculations. We focus on a weak, buoyant sphere in a stress weakening fluid subjected to far-field horizontal simple shear. The ascent velocity depends on two stress ratios: (1) the ratio of buoyancy stress to characteristic stress, controlling the transition from linear to power-law viscous flow, and (2) the ratio of regional stress associated with far-field shearing to characteristic stress. Comparing analytical estimates with numerical calculations, we find analytical estimates are accurate within a factor of two. However, discrepancies arise due to the analytical assumption that deviatoric stresses around the diapir are comparable to buoyancy stresses. Numerical results reveal significantly smaller deviatoric stresses. As deviatoric stresses govern stress-dependent, power-law, viscosity analytical estimates tend to overestimate stress weakening. We introduce a shape factor to improve accuracy. Additionally, we determine characteristic stresses for representative mantle and lower crustal flow laws and discuss practical implications in natural diapirism, such as sediment diapirs in subduction zones, magmatic plutons or exhumation of ultra-high-pressure rocks. Our study enhances understanding of diapir ascent velocities and associated stress conditions, contributing to a thorough comprehension of diapiric processes in geology.

Subsurface fluid injection stimulates complex hydromechanical interaction, necessitating the integration of geomechanical data across spatial and temporal scales to consider the sophisticated behavior. The induced seismic response is usually associated with reservoir architecture and pre-existing features that are three-dimensional, such as local stratigraphy, fractures, faults, and discontinuities. This study encompasses laboratory characterization of coupled hydromechanical response of rock cores extracted from reservoir - Mt. Simon sandstone, basal seal - Argenta sandstone, and crystalline basement - Precambrian rhyolite - formations in the Illinois Basin. High-resolution numerical modeling allows to consider the three-dimensional complexity of the injection site for Illinois Basin Decatur Project with spatial resolution comparable to one of the active seismic surveys. The laboratory-based porosity-permeability relation is combined with a three-dimensional porosity distribution developed from an inversion of active seismic data resulting in a detailed reconstruction of the evolution of the state of stress in formations where stress measurements are not performed. It appears that the microseismic clusters, mainly observed in the crystalline basement during the injection, are linked to zones experiencing more critically stressed conditions prior to injection. These zones have a potential for reactivation during the injection and are attributed to the specific local stratigraphy of the injection site, as well as transfer of triggering perturbations during the injection.

Deformation at tectonic plate boundaries involves coupling between rock deformation, fluid flow and metamorphic reactions, but quantifying this coupling is still elusive. We present a new two-dimensional hydro-mechanical-chemical numerical model and investigate the coupling between heterogeneous rock deformation and metamorphic (de)hydration reactions. Rock deformation consists of linear viscous compressible and power-law viscous shear deformation. Fluid flow follows Darcys law with a Kozeny-Carman type permeability. We consider a closed isothermal system and the reversible (de)hydration reaction: periclase and water yields brucite. In the models, fluid pressure within a circular or elliptical inclusion is initially below the periclase-brucite reaction pressure, and above in the surrounding. Inclusions exhibit a shear viscosity thousand times smaller than for the surrounding, because we assume that periclase-water and brucite regions have different effective viscosities. In models with circular inclusions, solid deformation has a minor impact on the evolution of fluid pressure, porosity and reaction front. Models with elliptical inclusions and far-field shortening generate higher rock pressure inside the inclusion compared to circular inclusions, and show a faster reaction-front propagation. The propagating reaction-front increases the inclusion surface and causes an effective, reaction-induced weakening of the heterogeneous rock. Weakening evolves strongly nonlinear with progressive strain. Distributions of fluid and rock pressure as well as magnitudes and directions of fluid and solid velocities are significantly different. The models mimic basic features of shear zones and plate boundaries and suggest a strong impact of heterogeneous rock deformation on (de)hydration reactions and associated reaction-induced weakening. The applied MATLAB algorithm is provided.

Melt transport across the ductile mantle is essential for oceanic crust formation or intraplate volcanism. Metasomatic enrichment of the lithospheric mantle demonstrates that melts chemically interact with the lithosphere. However, mechanisms of melt migration and the coupling of physical and chemical processes remain unclear. Here, we present a new thermo-hydro-mechanical-chemical (THMC) model for melt migration coupled to chemical differentiation. We study melt migration by porosity waves and consider a simple chemical system of forsterite-fayalite-silica. We solve the one-dimensional (1D) THMC model numerically using the finite difference method. Variables, such as solid and melt densities or and mass concentrations, are fully variable and functions of pressure (P), temperature (T) and total silica mass fraction (CTSiO2). These variables are pre-computed with thermodynamic Gibbs energy minimisation, which shows that dependencies of these variables to variations in P, T and CTSiO2 are considerably different. These P-T-CTSiO2 dependencies are implemented in the THMC model via parameterized equations. We consider P and T conditions relevant around the lithosphere-asthenosphere boundary and employ adiabatic and conductive geotherms. Variation of CTSiO2 changes the densities of solid and melt and has a strong impact on melt migration. We perform systematic 1D simulations to quantify the impact of initial distributions of porosity and CTSiO2 on the melt velocity. An adiabatic gradient generates higher melt velocities. Reasonable values for porosity, permeability, melt and compaction viscosities provide melt velocities between 10 [cm·yr-1] and 100 [m·yr-1]. We further discuss preliminary results of two 2D simulations showing blob-like and channel-like porosity waves.

Biot's equations describe the physics of hydro-mechanically coupled systems establishing the widely recognized theory of poroelasticity. This theory has a broad range of applications in Earth and biological sciences as well as in engineering. The numerical solution of Biot's equations is challenging because wave propagation and fluid pressure diffusion processes occur simultaneously but feature very different characteristic time scales. Analogous to geophysical data acquisition, high resolution and three dimensional numerical experiments lately re-defined state of the art. Tackling high spatial and temporal resolution requires a high-performance computing approach. We developed a multi-GPU numerical application to resolve the anisotropic elastodynamic Biot's equations that relies on a conservative numerical scheme to simulate, in a few seconds, wave fields for spatial domains involving more than 1.5 billion grid cells. We present a comprehensive dimensional analysis reducing the number of material parameters needed for the numerical experiments from ten to four. Furthermore, the dimensional analysis emphasizes the key material parameters governing the physics of wave propagation in poroelastic media. We perform a dispersion analysis as function of dimensionless parameters leading to simple and transparent dispersion relations. We then benchmark our numerical solution against an analytical plane wave solution. Finally, we present several numerical modeling experiments, including a three-dimensional simulation of fluid injection into a poroelastic medium. We provide the Matlab, symbolic Maple, and GPU CUDA C routines to reproduce the main presented results. The high efficiency of our numerical implementation makes it readily usable to investigate three-dimensional and high-resolution scenarios of practical applications.

We developed a numerical thermodynamics laboratory called “Thermolab” to study the effects of the thermodynamic behavior of non-ideal solution models on reactive transport processes in open systems. The equations of state of internally consistent thermodynamic datasets are implemented in MATLAB functions and form the basis for calculating Gibbs energy. A linear algebraic approach is used in Thermolab to compute Gibbs energy of mixing for multi-component phases to study the impact of the non-ideality of solution models on transport processes. The Gibbs energies are benchmarked with experimental data, phase diagrams and other thermodynamic software. Constrained Gibbs minimization is exemplified with MATLAB codes and iterative refinement of composition of mixtures may be used to increase precision and accuracy. All needed transport variables such as densities, phase compositions, and chemical potentials are obtained from Gibbs energy of the stable phases after the minimization in Thermolab. We demonstrate the use of precomputed local equilibrium data obtained with Thermolab in reactive transport models. In reactive fluid flow the shape and the velocity of the reaction front vary depending on the non-linearity of the partitioning of a component in fluid and solid. We argue that non-ideality of solution models has to be taken into account and further explored in reactive transport models. Thermolab Gibbs energies can be used in Cahn-Hilliard models for non-linear diffusion and phase growth. This presents a transient process towards equilibrium and avoids computational problems arising during precomputing of equilibrium data.