A two‐dimensional numerical simulation of lithospheric shortening shows the formation of a stable crustal‐scale shear zone due to viscous heating. The shear zone thickness is controlled by ...thermomechanical coupling that is resolved numerically inside the shear zone. Away from the shear zone, lithospheric deformation is dominated by pure shear, and tectonic overpressure (i.e., pressure larger than the lithostatic pressure) is proportional to the deviatoric stress. Inside the shear zone, deformation is dominated by simple shear, and the deviatoric stress decreases due to thermal weakening of the viscosity. To maintain a constant horizontal total stress across the weak shear zone (i.e., horizontal force balance), the pressure in the shear zone increases to compensate the decrease of the deviatoric stress. Tectonic overpressure in the weak shear zone can be significantly larger than the deviatoric stress at the same location. Implications for the geodynamic history of tectonic nappes including high‐pressure/ultrahigh‐pressure rocks are discussed.
Key Points
Viscous heating generates weak and stable crustal shear zones during shorteningHigh tectonic overpressure occurs in shear zone despite low differential stressHP‐UHP rocks can form in weak crustal thrust‐type shear zones
Large volumes of greenhouse gases such as CH
4 and CO
2 form by contact metamorphism of organic-rich sediments in aureoles around sill intrusions in sedimentary basins. Thermogenic gas generation and ...dehydration reactions in shale are treated numerically in order to quantify basin-scale devolatilization. We show that aureole thicknesses, defined as the zone of elevated metamorphism relative to the background level, vary within 30–250% of the sill thickness, depending on the temperature of the host-rock and intrusion, besides the sill thickness. In shales with total organic carbon content of >5
wt.%, CH
4 is the dominant volatile (85–135
kg/m
3) generated through organic cracking, relative to H
2O-generation from dehydration reactions (30–110
kg/m
3). Even using conservative estimates of melt volumes, extrapolation of our results to the scale of sill complexes in a sedimentary basin indicates that devolatilization can have generated ∼2700–16200
Gt CH
4 in the Karoo Basin (South Africa), and ∼600–3500
Gt CH
4 in the Vøring and Møre basins (offshore Norway). The generation of volatiles is occurring on a time-scale of 10–1000
years within an aureole of a single sill, which makes the rate of sill emplacement the time-constraining factor on a basin-scale. This study demonstrates that thousands of gigatons of potent greenhouse gases like methane can be generated during emplacement of Large Igneous Provinces in sedimentary basins.
Melt transport across the ductile mantle is essential for oceanic crust formation or intraplate volcanism. However, mechanisms of melt migration and associated chemical interaction between melt and ...solid mantle remain unclear. Here, we present a thermo‐hydro‐mechanical‐chemical (THMC) model for melt migration coupled to chemical differentiation. We consider melt migration by porosity waves and a 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 MgO and SiO2 mass concentrations, are functions of pressure (P), temperature (T), and total silica mass fraction (CTSiO2 ${C}_{\mathrm{T}}^{{\text{SiO}}_{2}}$). These variables are pre‐computed with Gibbs energy minimization and their variations with evolving P, T, and CTSiO2 ${C}_{\mathrm{T}}^{{\text{SiO}}_{2}}$ are implemented in the THMC model. We consider P and T conditions relevant around the lithosphere‐asthenosphere boundary. Systematic 1D simulations quantify the impact of initial distributions of porosity and CTSiO2 ${C}_{\mathrm{T}}^{{\text{SiO}}_{2}}$ on the melt velocity. Larger perturbations of CTSiO2 ${C}_{\mathrm{T}}^{{\text{SiO}}_{2}}$ cause larger melt velocities. An adiabatic or conductive geotherm cause fundamentally different vertical variations of densities and concentrations, and an adiabatic geotherm generates higher melt velocities. We quantify differences between melt transport (considering incompatible tracers), major element transport and porosity evolution. Melt transport is significant in the models. We also quantify the relative importance of four porosity variation mechanisms: (a) mechanical compaction and decompaction, (b) density variation, (c) compositional variation, and (d) solid‐melt mass exchange. In the models, (de)compaction dominates the porosity variation. We further discuss preliminary results of 2D THMC simulations showing blob‐like and channel‐like porosity waves.
Plain Language Summary
Melt transport across the ductile mantle is essential for oceanic crust formation or intraplate volcanism. However, melt transport processes are still incompletely understood and poorly quantified. Xenoliths (inclusions in igneous rock entrained during magma ascent) sampled by kimberlites (an igneous rock that erupted from the mantle) or intraplate basalts provide evidence that there is a reaction (metasomatism) between the rising melt and the solid mantle. However, the impact of chemical processes on melt migration remains unclear. Here, we present a new mathematical two‐phase (fluid‐solid) model, based on fundamental laws of physics and thermodynamics, which couples melt migration with chemical processes. We study melt migration around the lithosphere‐asthenosphere boundary and consider the solid (not molten) rocks as highly viscous fluids due to the high temperatures in these regions (at 80–100 km depth). We present one‐dimensional results of computer simulations and show that the variation of chemical components, such as silicon dioxide, changes the densities of the solid and melt, and can, hence, have a considerable impact on melt migration. We also present two‐dimensional simulations, which show the channelization of the rising melt.
Key Points
We present a thermo‐hydro‐mechanical‐chemical model for melt migration by porosity waves coupled to thermodynamic P–T–X phase calculations
Adiabatic or conductive geotherms cause different vertical variations of densities and concentrations as well as different melt velocities
Variations in Si concentration change solid and melt densities, trigger porosity waves and impact melt velocity, porosity and permeability
Common extrusion‐type models for the high‐ to ultrahigh‐pressure Adula nappe require a major normal fault along the top of this unit which is not conveyed in the structural record. This implies that ...such a normal fault existed but was completely erased during later deformational stages. However, there is evidence that decompression occurred during top‐to‐the‐foreland thrusting. We performed a new, purely structural kinematic restoration of the central part of the NFP20‐East cross section in order to estimate the burial depths of individual units without converting petrological pressure data into depth under the assumption that pressures were lithostatic. The results show that pressures within most of the units were close to but somewhat higher than lithostatic for several stages of the tectono‐metamorphic history. Only for the maximum burial stage of the Adula nappe, we estimate local tectonic overpressures of 40 to 80% of the lithostatic pressures. Accepting such an amount of overpressure, which is moderate compared to values theoretically possible, the Adula nappe was probably not subducted to subcrustal depth. We propose that the structural record of the Penninic nappe stack is quite complete and suggest that the decay of tectonic overpressure is a feasible explanation for decompression from eclogite‐ to amphibolite‐facies conditions during thrusting. Consequently, exhumation and convergence rates of the Eocene to Oligocene Alps may be smaller than previously assumed.
Key Points
Kinematics of Eastern Penninic Alps is reconstructed by structural record only
Tectonic overpressure solves conflict between structural and metamorphic record
Local and transient overpressure may determine P‐T‐t record of Adula (U)HP 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. We consider linear viscous compressible and power‐law viscous shear deformation. Fluid flow follows Darcy's law with a Kozeny‐Carman type permeability. We consider a closed isothermal system and the reversible (de)hydration reaction: periclase and water yields brucite. 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. For circular inclusions, solid deformation has a minor impact on the evolution of fluid pressure, porosity, and reaction front. For elliptical inclusions and far‐field shortening, rock pressure inside the inclusion is higher compared to circular inclusions, and reaction‐front propagation is faster. The reaction generates a large change in porosity (∼0.1%–55%) and in solid density (∼2,300–3,500 kg m−3), and the reaction front exhibits steep gradients of fluid pressure and porosity. Reaction‐front propagating increases the weak inclusion's surface and causes an effective, reaction‐induced weakening of the heterogeneous rock. Weakening evolves 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 show the importance of coupling deformation and metamorphism. The applied MATLAB algorithm is provided.
Plain Language Summary
Geodynamic processes at tectonic plate boundaries are complicated because rock deformation, fluid flow, and chemical reactions occur simultaneously. Investigating these coupled processes by direct observations is usually not possible, and investigating them with laboratory experiments is often not feasible. Alternatively, these coupled processes can be investigated with computer simulations. Here, we present a new two‐dimensional hydro‐mechanical‐chemical computer model to investigate the coupling of these processes. We consider a simple and reversible (de)hydration reaction: periclase (magnesium oxide) and water yields brucite (magnesium hydroxide). The initial fluid pressure within a circular or elliptical inclusion is initially below the periclase‐brucite reaction pressure, while in the surrounding it is above. Inclusions in the deforming rock are mechanically weaker than the surrounding. Models with elliptical inclusions generate higher rock pressure inside the inclusion compared to circular inclusions and show a faster reaction‐front propagation. The propagating reaction‐front causes an effective, reaction‐induced weakening of the heterogeneous rock. 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 suggest a strong impact of heterogeneous rock deformation on (de)hydration reactions and associated rheological weakening. The applied MATLAB algorithm is provided.
Key Points
Two‐dimensional (2D) hydro‐mechanical‐chemical model couples rock deformation, porous fluid flow, and reactions for brucite‐periclase (de)hydration reaction
Rock heterogeneities affect reaction‐front propagation and reaction‐induced weakening and cause differences in fluid and rock pressure
MATLAB code for numerical 2D hydro‐mechanical‐chemical model is provided
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: (a) the ratio of buoyancy stress to characteristic stress, controlling the transition from linear to power‐law viscous flow, and (b) 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.
Plain Language Summary
A diapir is a volume of rock that rises within a larger, denser rock mass due to its lower density and the force of gravity. Understanding the speed at which diapirs ascend is crucial for determining their significance in specific geologic settings, such as subduction zones. In this study, we use advanced computer simulations to calculate the ascent velocity of a spherical diapir within a denser surrounding material. The surrounding material is subjected to horizontal shearing, and its behavior resembles that of a nonlinear fluid, where its resistance to shear, known as viscosity, depends on the applied stress. By conducting three‐dimensional computer simulations, we not only test the accuracy of existing mathematical equations commonly used to estimate diapir velocity but also make improvements to enhance their precision. These equations help us estimate how quickly diapirs rise in different geodynamic environments. By advancing our understanding of diapir ascent velocities, we gain valuable insights into the processes that shape our planet's geological features.
Key Points
3D GPU‐based numerical calculations of diapir velocities in power‐law viscous fluid under far‐field stress
New analytical velocity estimates are controlled by two stress ratios and agree with numerical results
Stress weakening in tectonically active regions can increase diapir velocity by several orders of magnitudes
We present analytical derivations and 2-D numerical simulations that quantify magnitudes of deviatoric stress and tectonic overpressure (i.e. difference between the pressure, or mean stress, and the ...lithostatic pressure) by relating them to lateral variations in the gravitational potential energy (GPE). These predictions of tectonic overpressure and deviatoric stress associated with GPE differences are independent of rock rheology (e.g. viscous or elastic) and rock strength. We consider a simple situation with lowlands and mountains (plateau). We use a numerical two-layer model consisting of a crust with higher Newtonian viscosity than that in the mantle, and also a three-layer model in which the two-layer lithosphere overlies a much less viscous asthenosphere. Our results (1) explain why estimates for the magnitude of stresses in Tibet, previously published by different authors, vary by a factor of two, (2) are applied to test the validity of the thin sheet approximation, (3) show that the magnitude of the depth-integrated tectonic overpressure is equal to the magnitude of the depth-integrated deviatoric stress if depth-integrated shear stresses on vertical and horizontal planes within the lithosphere are negligible (the thin sheet approximation) and (4) show that under thin sheet approximation tectonic overpressure is required to build and support continental plateaus, such as in Tibet or in the Andes, even if the topography and the crustal root are in isostatic equilibrium. Under thin sheet approximation, the magnitude of the depth-integrated tectonic overpressure is equal to the depth-integrated horizontal deviatoric stress, and both are approximately 3.5 × 1012 N m−1 for Tibet. The horizontal driving force per unit length related to lateral GPE variations around Tibet is composed of the sum of both tectonic overpressure and deviatoric stress, and is approximately 7 × 1012 N m−1. This magnitude exceeds previously published estimates for the force per unit length required to fold the Indo-Australian Plate south of India, and hence the uplift of the Tibetan plateau could have folded the Indian Plate. We also discuss the mechanical conditions that are necessary to achieve isostasy, for which the lithostatic pressure is constant at a certain depth. The results show that tectonic overpressure can exist at a certain depth even if all deviatoric stresses are zero at this depth, because this tectonic overpressure is related to horizontal gradients of vertical shear stresses integrated across the entire depth of the lithosphere. The magnitude of the depth-integrated tectonic overpressure of 3.5 × 1012 N m−1 implies that the pressure estimated from observed mineral assemblages in crustal rocks is likely significantly different from the lithostatic pressure, and pressure recorded by crustal rocks is not directly related to depth. In case of significant weakening of the entire lithosphere by any mechanism our analytical and numerical studies provide a simple estimation of tectonic overpressure via variations in GPE.
Focused fluid flow is common in sedimentary basins worldwide, where flow structures often penetrate through sandy reservoir rocks, and clay‐rich caprocks. To better understand the mechanisms forming ...such structures, the impacts of the viscoelastic deformation and strongly nonlinear porosity‐dependent permeability of clay‐rich materials are assessed from an experimental and numerical modeling perspective. The experimental methods to measure the poroviscoelastic and transport properties of intact and remolded shale have been developed, and the experimental data is used to constrain the numerical simulations. It is demonstrated that viscoelastic deformation combined with nonlinear porosity‐dependent permeability triggers the development of localized flow channels, often imaged as seismic chimneys. The permeability inside a channel increases by several orders of magnitude compared to the background values. In addition, the propagation time scale and the channel size strongly depend on the material properties of the fluid and the rock. The time‐dependent behavior of the clay‐rich rock may play a key role in the long‐term integrity of the subsurface formations.
Key Points
Clay‐rich rocks exhibit bulk viscous deformation in laboratory experiments
Focused fluid flow might develop through clay‐rich rock with implications for injectivity and leakage risk
Model simulations based on experimental parameters give estimates for expected leakage rates
We present two-dimensional numerical simulations of lithospheric shortening with a crust containing weak and strong inclusions. Thermo-mechanical coupling is included, and a crustal-scale shear zone ...develops self-consistently due to viscous heating and thermal softening of temperature dependent viscosities. Several tests for crustal conditions are performed showing that 1) the thickness of and strain rates within the shear zone are independent on the numerical resolution and applied numerical method (finite element and finite difference method), 2) the shear zone is stable and rotates during large strain deformation, 3) the numerical algorithm conserves total thermal and mechanical energies, and 4) the bulk horizontal force balance is fulfilled during large strain deformation. A fold nappe develops around the shear zone in the lithospheric shortening simulation. In this simulation the stresses in the crust are limited by a friction angle of 30°. Significant tectonic overpressure (PO) occurs in strong lower crustal rocks and in strong inclusions. Significant PO also occurs in a weak inclusion that is only partly surrounded by strong crustal rock suggesting that a continuous strong “vessel” is not required to generate significant PO in weak rocks. Maximal values of PO are ~2.2GPa with corresponding deviatoric stresses ~1.5GPa and occur in a depth of ~42km. Maximal pressure of~3.4GPa and maximal temperatures>700°C occur during the formation of the fold nappe in crustal depth. Synthetic pressure–temperature paths exhibit entire cycles of pressure and temperature increase and decrease, and suggest that crustal rocks in depths<50km can reach the ultrahigh pressure metamorphic facies fields. Applications to tectonic nappes with high and ultra-high pressure rocks in the Western Alps are discussed, and a dynamic model for the evolution of fold nappes in the Western Alps is proposed.
•A new dynamic model for tectonic nappes in the Western Alps is presented.•High and ultra-high pressures can occur in crustal depth due to tectonic overpressure.•Thickness of shear zones due to shear heating is independent on numerical resolution.
PDF
