The exponential flow law applied to necking and folding of a ductile layer

The uniaxial exponential law, , has been applied to experimental results for steady-state creep, including that of wet and dry olivine, pyroxenite and carbonates, at a deviatoric stress greater than ∼100 MPa. Such stress levels likely occur in the upper-mantle lithosphere and middle crust. In a layered rock, in layer-parallel extension or shortening, high deviatoric stress can occur in the stiffest layers, for example, a dolomite layer in fine-grained marble. Under assumptions of isotropy and incompressibility, the uniaxial law yields where J2 = sijsij/2 is an invariant of the deviatoric stress tensor, sij. Linearization about a homogeneous basic-state of flow yields equations identical in form to those obtained for the familiar power law. This establishes expressions for an effective viscosity, and stress exponent, , where is an invariant of the basic-state deformation rate, . These results allow application of existing analytical folding and necking solutions for a rock layer of power-law fluid to a rock layer with an exponential flow law. The effective stress exponent for the exponential flow law increases with decreasing temperature, through the dependence of C on the latter and, weakly, with increasing deformation rate. For dry olivine, effective stress exponents are between 10 and 30 for temperatures between 400 and 600°C with little dependence on deformation rate. Finite element simulations employing full non-linear forms of the flow laws show that large strain necking is nearly identical for power law and exponential flow laws. The results suggest that the instability in necking and folding in ductile rock layers can be considerably stronger than inferred from results based on flow laws representing diffusion and dislocation creep. The large values of the effective stress exponent, ne > ∼15, that may be attained for exponential flow laws can account for observed outcrop-scale ductile neckin

Spectral modification of seismic waves propagating through solids exhibiting a resonance frequency: a 1-D coupled wave propagation-oscillation model

A 1-D model is presented that couples the microscale oscillations of non-wetting fluid blobs in a partially saturated poroelastic medium with the macroscale wave propagation through the elastic skeleton. The fluid oscillations are caused by surface tension forces that act as the restoring forces driving the oscillations. The oscillations are described mathematically with the equation for a linear oscillator and the wave propagation is described with the 1-D elastic wave equation. Coupling is done using Hamilton's variational principle for continuous systems. The resulting linear system of two partial differential equations is solved numerically with explicit finite differences. Numerical simulations are used to analyse the effect of solids exhibiting internal oscillations, and consequently a resonance frequency, on seismic waves propagating through such media. The phase velocity dispersion relation shows a higher phase velocity in the high-frequency limit and a lower phase velocity in the low-frequency limit. At the resonance frequency a singularity in the dispersion relation occurs. Seismic waves can initiate oscillations of the fluid by transferring energy from solid to fluid at the resonance frequency. Due to this transfer, the spectral amplitude of the solid particle velocity decreases at the resonance frequency. After initiation, the oscillatory movement of the fluid continuously transfers energy at the resonance frequency back to the solid. Therefore, the spectral amplitude of the solid particle velocity is increased at the resonance frequency. Once initiated, fluid oscillations decrease in amplitude with increasing time. Consequently, the spectral peak of the solid particle velocity at the resonance frequency decreases with tim

Impact of upper mantle convection on lithosphere hyperextension and subsequent horizontally forced subduction initiation

Many plate tectonic processes, such as subduction initiation, are embedded in long-term (> 100 Myr) geodynamic cycles often involving subsequent phases of extension, cooling without plate deformation and convergence. However, the impact of upper mantle convection on lithosphere dynamics during such long-term cycles is still poorly understood. We have designed two-dimensional upper mantle- scale (down to a depth of 660 km) thermo-mechanical numerical models of coupled lithosphere–mantle deformation. We consider visco–elasto–plastic deformation including a combination of diffusion, dislocation and Peierls creep law mechanisms. Mantle densities are calculated from petrological phase diagrams (Perple_X) for a Hawaiian pyrolite. Our models exhibit realistic Rayleigh numbers between 1e6 and 1e7, and the model temperature, density and viscosity structures agree with geological and geophysical data and observations. We tested the impact of the viscosity structure in the asthenosphere on upper mantle convection and lithosphere dynamics. We also compare models in which mantle convection is explicitly modelled with models in which convection is parameterized by Nusselt number scaling of the mantle thermal conductivity. Further, we quantified the plate driving forces necessary for subduction initiation in 2D thermo- mechanical models of coupled lithosphere–mantle deformation. Our model generates a 120 Myr long geodynamic cycle of subsequent extension (30 Myr), cooling (70 Myr) and convergence (20 Myr) coupled to upper mantle convection in a single and continuous simulation. Fundamental features such as the formation of hyperextended margins, upper mantle convective flow and subduction initiation are captured by the simulations presented here. Compared to a strong asthenosphere, a weak asthenosphere leads to the following differences: smaller value of plate driving forces necessary for subduction initiation (15 TN/m instead of 22 TN/m) and locally larger suction forces. The latter assists in establishing single-slab subduction rather than double-slab subduction. Subduction initiation is horizontally forced, occurs at the transition from the exhumed mantle to the hyperextended passive margin and is caused by thermal softening. Spontaneous subduction initiation due to negative buoyancy of the 400 km wide, cooled, exhumed mantle is not observed after 100 Myr in model history. Our models indicate that long-term lithosphere dynamics can be strongly impacted by sub-lithosphere dynamics. The first-order processes in the simulated geodynamic cycle are applicable to orogenies that resulted from the opening and closure of embryonic oceans bounded by magma-poor hyperextended rifted margins, which might have been the case for the Alpine orogeny

Numerical Simulations Reproduce Field Observations Showing Transient Weakening During Shear Zone Formation by Diffusional Hydrogen Influx and H2O Inflow

Exposures on Holsnøy island (Bergen Arcs, Norway) indicate fluid infiltration through fractures into a dry, metastable granulite, which triggered a kinetically delayed eclogitization, a transient weakening during fluid-rock interaction, and formation of shear zones that widened during shearing. It remains unclear whether the effects of grain boundary-assisted aqueous fluid inflow on the duration of granulite hydration were influenced by a diffusional hydrogen influx accompanying the fluid inflow. To better estimate the fluid infiltration efficiencies and the parameter interdependencies, a 1D numerical model of a viscous shear zone is utilized and validated using measured mineral phase abundance distributions and H2O-contents in nominally anhydrous minerals of the original granulite assemblage to constrain the hydration by aqueous fluid inflow and diffusional hydrogen influx, respectively. Both hydrations are described with a diffusion equation and affect the effective viscosity. Shear zone kinematics are constrained by the observed shear strain and thickness. The model fits the phase abundance and H2O-content profiles if the effective hydrogen diffusivity is approximately one order of magnitude higher than the diffusivity for aqueous fluid inflow. The observed shear zone thickness is reproduced if the viscosity ratio between dry granulite and deforming, reequilibrating eclogite is ∼104 and that between dry granulite and hydrated granulite is ∼102. The results suggest shear velocities <10−2 cm/a, hydrogen diffusivities of ∼10−13±1 m2/s, and a shearing duration of <10 years. This study successfully links and validates field data to a shear zone model and highlights the importance of hydrogen diffusion for shear zone widening and eclogitization

A simple computer program for calculating stress and strain rate in 2D viscous inclusion-matrix systems

Computer-based numerical solutions of geomechanical problems are important to understand the processes forming rock structures as well as to quantify the associated pressure, stresses and strain rates. However, the development of such computer programs and the underlying mathematical methods are commonly not taught in a standard structural geology curriculum. Here, we present a simple computer program to calculate the stress, pressure, velocity and strain rate fields for two-dimensional (2D) viscous inclusion-matrix systems under pure shear and simple shear. The main aim of our contribution is to explain this computer program in a simple and transparent way, so that it can be used for introductory courses on geomechanical numerical modelling in structural geology. We present the governing equations of 2D viscous deformation and program the equations in the same order and style, so that the equations are still recognizable in the computer program. The computer program can treat stiff and weak inclusions of any shape and considers both linear and power-law viscous flow laws. We present numerical calculations for various inclusion-matrix scenarios. The program is written with the software MATLAB, is provided as supplementary material, and can also be run with the freely available software GNU Octave

Relationship between tectonic overpressure, deviatoric stress, driving force, isostasy and gravitational potential energy

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 GP

A spectral/finite difference method for simulating large deformations of heterogeneous, viscoelastic materials

A numerical algorithm is presented that simulates large deformations of heterogeneous, viscoelastic materials in two dimensions. The algorithm is based on a spectral/finite difference method and uses the Eulerian formulation including objective derivatives ofthe stress tensor in the rheological equations. The viscoelastic rheology is described bythe linear Maxwell model, which consists of an elastic and viscous element connected inseries. The algorithm is especially suitable to simulate periodic instabilities. The derivatives in the direction of periodicity are approximated by spectral expansions, whereas the derivatives in the direction orthogonal to the periodicity are approximated by finite differences. The 1‐D Eulerian finite difference grid consists of centre and nodal points and has variable grid spacing. Time derivatives are approximated with finite differences using an implicit strategy with a variable time step. The performance of the numerical code is demonstrated by calculation, for the first time, of the pressure field evolution during folding of viscoelastic multilayers. The algorithm is stable for viscosity contrasts up to 5 × 105, which demonstrates that spectral methods can be used to simulate dynamical systems involving large material heterogeneities. The successful simulations show that combined spectral/finite difference methods using the Eulerian formulation are a promising tool to simulate mechanical processes that involve large deformations, viscoelastic rheologies and strong material heterogeneitie

Melt Migration and Chemical Differentiation by Reactive Porosity Waves

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 (urn:x-wiley:15252027:media:ggge22741:ggge22741-math-0001). These variables are pre-computed with Gibbs energy minimization and their variations with evolving P, T, and urn:x-wiley:15252027:media:ggge22741:ggge22741-math-0002 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 urn:x-wiley:15252027:media:ggge22741:ggge22741-math-0003 on the melt velocity. Larger perturbations of urn:x-wiley:15252027:media:ggge22741:ggge22741-math-0004 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

Buoyancy versus shear forces in building orogenic wedges

International audienceThe dynamics of growing collisional orogens are mainly controlled by buoyancy and shear forces. However, the relative importance of these forces, their temporal evolution and their impact on the tectonic style of orogenic wedges remain elusive. Here, we quantify buoyancy and shear forces during collisional orogeny and investigate their impact on orogenic wedge formation and exhumation of crustal rocks. We leverage two-dimensional petrological–thermomechanical numerical simulations of a long-term (ca. 170 Myr) lithosphere deformation cycle involving subsequent hyperextension, cooling, convergence, subduction and collision. Hyperextension generates a basin with exhumed continental mantle bounded by asymmetric passive margins. Before convergence, we replace the top few kilometres of the exhumed mantle with serpentinite to investigate its role during subduction and collision.We study the impact of three parameters: (1) shear resistance, or strength, of serpentinites, controlling the strength of the evolving subduction interface; (2) strength of the continental upper crust; and (3) density structure of the subducted material. Densities are determined by linearized equations of state or by petrological-phase equilibria calculations. The three parameters control the evolution of the ratio of upward-directed buoyancy force to horizontal driving force, FB/FD=ArF, which controls the mode of orogenic wedge formation: ArF≈0.5 causes thrust-sheet-dominated wedges, ArF≈0.75 causes minor wedge formation due to relamination of subducted crust below the upper plate, and ArF≈1 causes buoyancy-flow- or diapir-dominated wedges involving exhumation of crustal material from great depth (>80 km). Furthermore, employing phase equilibria density models reduces the average topography of wedges by several kilometres.We suggest that during the formation of the Pyrenees ArF⪅0.5 due to the absence of high-grade metamorphic rocks, whereas for the Alps ArF≈1 during exhumation of high-grade rocks and ArF⪅0.5 during the post-collisional stage. In the models, FD increases during wedge growth and subduction and eventually reaches magnitudes (≈18 TN m−1) which are required to initiate subduction. Such an increase in the horizontal force, required to continue driving subduction, might have “choked” the subduction of the European plate below the Adriatic one between 35 and 25 Ma and could have caused the reorganization of plate motion and subduction initiation of the Adriatic plate