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.

Serpentinite subduction and associated dehydration vein formation are important for subduction zone dynamics and water cycling. Field observations suggest that en échelon olivine veins in serpentinite mylonites formed by dehydration during simultaneous shearing of serpentinite. Here, we test a hypothesis of shear-driven formation of dehydration veins with a two-dimensional hydro-mechanical-chemical numerical model. We consider the reaction antigorite + brucite = forsterite + water. Shearing is viscous and the shear viscosity decreases with increasing porosity. Total and fluid pressures are initially homogeneous and in the serpentinite stability field. Initial perturbations in porosity, and hence viscosity, cause fluid pressure perturbations during simple shearing. Dehydration nucleates where fluid pressure decreases locally below the thermodynamic pressure defining the reaction boundary. During shearing, dehydration veins grow in direction parallel to the maximum principal stress and serpentinite transforms into olivine inside the veins. Simulations show that the relation between compaction length and porosity as well as the ambient pressure have a strong impact on vein formation, while the orientation of the initial porosity perturbation and a pressure-insensitive yield stress have a minor impact. Porosity production associated with dehydration is controlled by three mechanisms: solid volumetric deformation, solid density variation and reactive mass transfer. Vein formation is self-limiting and slows down due to fluid flow decreasing fluid pressure gradients. We discuss applications to natural olivine veins as well as implications for slow slip and tremor, transient weakening, anisotropy generation and the formation of shear-driven high-porosity bands in the absence of a dehydration reaction.

Serpentinite subduction and the associated formation of dehydration veins is important for subduction zone dynamics and water cycling. Field observations suggest that en-échelon olivine veins in serpentinite mylonites formed by dehydration during simultaneous shearing of ductile serpentinite. Here, we test a hypothesis of shear-driven formation of dehydration veins with a two-dimensional hydro-mechanical-chemical numerical model. We consider the reaction antigorite + brucite = forsterite + water. Shearing is viscous and the shear viscosity decreases exponentially with porosity. The total and fluid pressures are initially homogeneous and in the antigorite stability field. Initial perturbations in porosity, and hence viscosity, cause fluid pressure perturbations. Dehydration nucleates where the fluid pressure decreases locally below the thermodynamic pressure defining the reaction boundary. Dehydration veins grow during progressive simple-shearing in a direction parallel to the maximum principal stress, without involving fracturing. The porosity evolution associated with dehydration reactions is controlled to approximately equal parts by three mechanisms: volumetric deformation, solid density variation and reactive mass transfer. The temporal evolution of dehydration veins is controlled by three characteristic time scales for shearing, mineral-reaction kinetics and fluid-pressure diffusion. The modelled vein formation is self-limiting and slows down due to fluid flow decreasing fluid pressure gradients. Mineral-reaction kinetics must be significantly faster than fluid-pressure diffusion to generate forsterite during vein formation. The self-limiting feature can explain the natural observation of many, small olivine veins and the absence of few, large veins. We further discuss implications for transient weakening during metamorphism and episodic tremor and slow-slip in subduction zones.

Double seismic zones (DSZs), parallel planes of intermediate-depth earthquakes inside oceanic slabs, have been observed in a number of subduction zones and may be a ubiquitous feature of downgoing oceanic plates. Focal mechanism observations from DSZ earthquakes sample the intraslab stress field at two distinct depth levels within the downgoing lithosphere. A pattern of downdip compressive over downdip extensive events was early on interpreted to indicate an unbending-dominated intraslab stress field. In the present study, we show that the intraslab stress field in the depth range of DSZs is much more variable than previously thought. Compiling DSZ locations and mechanisms from literature, we observe that the “classical’ pattern of compressive over extensive events is only observed at about half of the DSZ locations around the globe. The occurrence of extensional mechanisms across both planes accounts for most other regions. To obtain an independent estimate of the bending state of slabs at intermediate depths, we compute (un)bending estimates from slab geometries taken from the slab2 compilation of slab surface depths. We find no clear global prevalence of slab unbending at intermediate depths, and the occurrence of DSZ seismicity does not appear to be limited to regions of slab unbending. Focal mechanism observations are frequently inconsistent with (un)bending estimates from slab geometries, which may imply that bending stresses are not always prevalent, and that other stress types such as in-plane tension due to slab pull or shallow compression due to friction along the plate interface may also play an important role.

The evolution and distribution of metamorphic rocks throughout the western European Alps is indicative of subduction-related metamorphism. The present-day distribution of metamorphic rocks in the Western Alps exhibits a regional trend, with an internal high-pressure domain and decreasing grade towards the foreland. However, the processes by which high-grade continental rocks are formed and exhumed, as well as the evolution of the metamorphic architecture remains unclear. Here, we present a two-dimensional petrological-thermomechanical model to investigate the evolution and distribution of metamorphic facies within an orogenic wedge formed by subduction and continental collision. The model simulates an entire geodynamic cycle of extension, with passive margin formation and mantle exhumation, followed by thermal equilibration without applied far-field deformation, convergence, with subduction initiation, basin closure and collision. After thermal equilibration, we consider ad-hoc the serpentinization of the exhumed mantle. Models developing a weak subduction interface, due to 6 km serpentinite thickness, display a laterally varying peak metamorphic facies distribution, with the highest grade rocks within the core of the orogeny, agreeing with distributions in the Western Alps. In contrast, models with a stronger subduction interface (3 km serpentinite thickness) develop an orogenic wedge with a vertical metamorphic gradient. The metamorphic distribution is calculated using the peak P and T values of 10’000 numerical markers during their modelled P-T trajectories. The models indicate, during overall convergence, local extensional tectonics between the exhuming material and overriding plate, whereby the upper-plate hanging-wall is unroofed, moving with a normal sense of shear relative to the exhuming high-pressure rocks.

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.