Search

Radial microcracks surrounding retrogressed SiO2 inclusions in UHP garnets are key microstructural observations allowing to constrain the mechanisms of exhumation of ultra-high-pressure (UHP) rocks. The major challenge lies in identifying whether the microstructures formed during their ascent from mantle depths, or as a consequence of transient variations in the tectonic regime. By combining petrographic observations, petrochronological data and numerical thermo-mechanical modelling, we show that radial cracks around SiO2 inclusions in ultrahigh-pressure garnets from Dora Maira are caused by ultrafast decompression during the early stage of exhumation (< 0.5 Ma). Decompression rates higher than 10-14 s-1 are, for the first time, inferred from natural microstructures independently of current petrochronological estimates1. We demonstrate that the SiO2 phase transition generates shear stresses sufficiently large to trigger plastic yielding, resulting in the generation and propagation of radial and bent shear bands, mimicking the fractures observed in UHP garnet. Our results question the traditional interpretation of the exhumation from great depth of ultrahigh-pressure tectonic. Instead, we propose that such ultrafast decompression rates are related to transient changes in the stress state of the buried continental lithosphere, favoring an exhumation mechanism involving nappe stacking. 1 Rubatto, D. & Hermann, J. Exhumation as fast as subduction? Geology 29, 3–6 (2001).

The processes that govern rock (trans)formation (deposition, deformation, segregation, metamorphism) can result in the development of layering and rock fabrics. Rocks can thus exhibit extrinsic or intrinsic anisotropy at various spatial scales. Anisotropy has important mechanical consequences, in particular, for strain localisation in the lithosphere. This effect is typically not included in geodynamic models. Mechanical anisotropy can be modelled by explicitly modelled by numerically resolving layers of different strengths. Due to the expensive computational cost, this approach is not suitable for large scale geodynamic models. The latter may rather benefit from an upscaling approach that involves anisotropic constitutive laws. To model the evolution of such material Mühhlaus, (2002) proposed the use of the director vector which corresponds to a single orientation that is changing throughout the process of deformation. We have implemented visco-elasto-plastic anisotropic constitutive laws and the director vector approach in the geodynamic simulation tool MDoodz7.0. Here we present the rheological implementation, we show some simple simulations involving anisotropic flow and discuss the potential role of anisotropy for large-scale geodynamic processes.

The mechanical characteristics of a shell, having a double curvature, are fundamentally different to the characteristics of a plate, having no curvature in its undeformed state. Geometrically, the Earth’s lithosphere is a shell rather than a plate. However, most geodynamic numerical models applied to study the deformation of the lithosphere do not consider this curvature. It is currently unclear whether the shell-type geometry of the lithosphere has a significant impact on lithosphere deformation on the scale of few 1000 kilometers. This study investigates the importance of considering lithospheric shells and compares numerical results of a shortening shell-type and plate-type lithosphere. We apply the two-dimensional state-of-the-art thermo-mechanical code MDoodz (Duretz et al. 2021). We consider a shortening lithosphere in an initially curved and in an initially rectangular geometry and calculate the spatio-temporal stress distribution inside the deforming lithosphere. We further present preliminary results on the effects and relative importance of various softening mechanism, leading to strain localization and subduction initiation, such as thermal softening, grain size reduction, or anisotropy generation due to fabric development. REFERENCESDuretz T., R. de Borst and P. Yamato (2021), Modeling Lithospheric Deformation Using a Compressible Visco-Elasto-Viscoplastic Rheology and the Effective Viscosity Approach, Geochemistry, Geophysics, Geosystems, Vol. 22 (8), e2021GC009675

Eclogitization constitutes one of the most emblematic transformations in continental subduction zones, where conversion of initially dry lower crustal rocks into eclogite facies rocks correlates with the occurrence of seismogenic events. This reaction is generally considered to occur at high pressure conditions during hydration of dry granulite. Several models using « ad hoc » diffusion equation exist to model this hydration process and the consequences of reaction-induced changes in terms of rheology and density. However, to our knowledge, there is no quantitative model allowing to physically explain how fluids propagate inside a dry rock (i.e. with no porosity at all) and how reaction-induced alteration front widens over time. In this study, we therefore propose a new fully coupled hydro-chemical model wherein a two-phase flow model is coupled with the eclogitization reaction. We use a mass conservative approach, solving total mass and solid mass equations, in a closed isothermal system. Fluid and solid densities are calculated with lookup tables from equilibrium thermodynamics. Our model shows that a fluid pressure pulse generates a pressure gradient that can be associated with the densification reaction when the pressure required for the eclogitization is reached. This reaction generates a large increase in porosity (0 to ~16%) and subsequent porous fluid flow inside the initially dry granulite. This process is then sustained as long as the fluid pulse is maintained, and ends shortly after the fluid pressure pulse stops. However, high pressure within the reacted area can persist for a long period of time. A parametric study allowing to constrain both the duration and the widening of the reaction area is proposed as well as an application to the emblematic case study of the eclogitized granulites of Holsnoy (Bergen Arcs, Norway).

Eclogitization reactions in mafic rocks involve large volume changes, porosity evolution and fluid transfer. They impact many important geological processes such as the localization of deformation and fluid channeling at intermediate depth in subduction zone. The study of exhumed eclogitic bodies in orogens shows that eclogitization of the oceanic crust is heterogeneous from both a structural and metamorphic point of view. For example, in the European Alps, the Allalin metagabbro shows high strain areas, consisting of hydrous metagabbros, fully equilibrated under eclogite-facies conditions during the Alpine orogeny. Conversely, large volumes of low strain, fluid-undersaturated gabbros remained largely unaffected by the high-pressure (HP) metamorphism, locally preserving igneous textures and even, occasionally, relics of their magmatic mineralogy. The intensity of deformation as well as the degree of eclogitization in the metagabbro have been shown to be directly related to the extent of pre-Alpine hydration during high-temperature hydrothermal alteration [1]. However, the influence of this degree of hydration on (1) reaction kinetics and/or (2) enhancing rheological contrasts leading to heterogeneous deformation patterns and metamorphic conditions is still debated.In order to address this issue, we propose a multidisciplinary study involving petrographic and microtextural observations combined with 2D thermo-mechanical numerical models allowing to discuss the role of pre-Alpine hydrothermal alteration on the development of HP metamorphic assemblages.We present petrographic and textural data from three different types of rocks from the Allalin metagabbros: i) undeformed and mostly untransformed metagabbros, with relics of igneous augite and plagioclase, ii) coronites, with olivine pseudomorphs showing different levels of hydration, rimmed by a garnet corona, and iii) eclogitized metagabbros, with olivine and plagioclase sites fully replaced by high-pressure assemblages.The role of protolith hydration on the observed range in metamorphic facies is then tested by using 2D thermo-mechanical models that allow to simulate the deformation of a strong and dry rock with several randomly oriented weak and hydrous zones. Our results show that the shearing of heterogeneous rock can lead to the formation of localized ductile shear zone within a matrix that remains relatively undeformed but where plastic deformation can occur. The associated P field is also highly heterogeneous, with P ranging from 1 to 3 GPa. The deformation patterns and P modelled may suggest that locally hydrated portions of the gabbro acted as rheological perturbations sufficiently efficient in producing the structural and metamorphic record now observed in the field. [1] Barnicoat, A. C. & Cartwright, I. (1997) Journal of Metamorphic Geology 15, 93–104

Since their discovery in 1928, deep earthquakes have been the subject of extensive research to unravel their nucleation and rupture mechanisms. Due to the elevated pressures and temperatures at depths below 50 km, brittle failure becomes less likely and ductile deformation is favored. To date, there is no consensus on the mechanisms resulting in deep earthquake generation. Three main mechanisms (dehydration embrittlement, transformational faulting and thermal runaway) have been proposed to cause deep earthquakes, but neither of them has been sufficiently quantified to yield a definite answer under which conditions they are active.Here, we explore the feasibility of the thermal runaway hypothesis using 1D and 2D thermo-mechanical models. In particular, we investigate the impact of grain size reduction in conjunction with shear heating to see whether grain size reduction and shear heating are competitive mechanisms (which would prevent thermal runaway) or whether they are collaborative. Our results show that the combination of both mechanisms facilitates thermal runaway and significantly reduces the stress required for the occurrence of thermal runaway. We then investigate whether this combined failure mechanism may explain the seismicity observed in regions of detaching lithosphere, such as the Hindu Kush and the Vrancea earthquake nests.

The Face-Centered Finite Volume (FCFV) is a newly developed discretisation technique that has been applied to a variety of engineering problems. This approach is based on the hybridisable discontinuous Galerkin formulation with constant degree approximations. The FCFV is particularly attractive approach since it meets numerous essential criteria for successful geodynamic modelling. It offers full geometric flexibility, natural free surface boundary condition, second order accuracy velocity-field solutions, no oscillatory pressure modes, relatively low computational cost and adequate treatment of jump conditions at material interfaces. Here we present the implementation of Poisson and Stokes solvers in the Julia computing language. Here we present the implementation of Poisson and Stokes solvers using the performant Julia language. We discuss several solving strategies including direct-iterative and iterative pseudo-transient approaches, the latter executing efficiently on Graphical Processing Units. We extend the original FCFV Stokes formulation to account for discontinuous viscosity case and discuss the implementation of complex visco-elasto-plastic rheologies.

The formation and maintenance of narrow, lithospheric shear zones and their importance in plate-tectonics remain one of the major problems in geodynamics. While the cause and consequence of strain localization and weakening within the lithosphere remain debated, it is clear that these processes play an essential role in lithospheric deformation across a wide range of spatio-temporal scales. Here, we analyze the efficiency of strain localization in a 2-D visco-elasto-plastic medium for a strain-dependent weakening and healing (SDWH) rheology using 2-D numerical, thermo-mechanical experiments with kinematic boundary conditions. Such a parameterized rheology successfully mimics more complex transient weakening and healing processes, akin to a grain-size sensitive composite (diffusion and dislocation creep) rheology. In addition, the SDWH rheology allows for memory of deformation. This enables self-consistent formation and reactivation of inherited weak zones within the lithosphere and sustains those weak zones over an extended period of time. We further analyze the resulting shear zone patterns and seek to answer the questions: What is the typical, effective intensity of strain localization? What are the dimensions of the resulting shear zones? Are such shear zones mesh-dependent in numerical models and, if so, can we exploit existing regularization approaches for the SDWH rheology?

Metamorphic reactions can lead to drastic changes in rocks mechanical properties. Indeed, during such transformations, the nucleation of new phases with different strength, grain size and/or density compared to the primary phases is enhanced, and transient processes due to the ongoing reaction are then activated.Eclogitization of lower crustal rocks during continental subduction constitutes an emblematic transformation illustrating these processes. In such tectonic context, it has been shown that eclogitization seems to be closely associated with the occurrence of seismogenic events. However, the mechanisms that trigger brittle failure in such high pressure environments remain highly debated. Indeed, whether the change in density or the change in rheology can lead to embrittlement is still enigmatic.By using 2D compressible mechanical numerical models we studied the impact of the strong negative volume change of the eclogitization reaction on the rocks rheological behaviour. We show that eclogitization-induced density change occurring out of equilibrium can, by itself, generates sufficient shear stress to fail the rocks at high-pressure conditions.Rupture initiation at depth in continental subduction zones could therefore be explained by volume changes, even without considering the modifications of the rheological properties induced by the transformation. Our results also indicate that the negative volume change associated with brittle failure can enhance the propagation of the eclogitization process by a runaway mechanism as long as the reaction is not limited by the lack of reactants.

Strain localization and associated softening mechanisms in a deforming lithosphere are important for subduction initiation or the generation of tectonic nappes during orogeny. Many strain localization and softening mechanisms have been proposed as being important during the viscous, creeping, deformation of the lithosphere, such as thermal softening, grain size reduction, reaction-induced softening or anisotropy development. However, which localization mechanism is the controlling one and under which deformation conditions is still contentious. In this contribution, we focus on strain localization in viscous material due to the generation of anisotropy, for example due to the development of a foliation. We numerically model the generation and evolution of anisotropy during two-dimensional viscous simple shear in order to quantify the impact of anisotropy development on strain localization and on the effective softening. We calculate the finite strain ellipse during viscous deformation. The aspect ratio of the finite strain ellipse serves as proxy for the magnitude and evolution of anisotropy, which determines the ratio of normal to tangential viscosity. To track the orientation of the anisotropy during deformation we apply a director method. We benchmark our implementation of anisotropy by comparing results of two independently developed numerical algorithms based on the finite difference method: one algorithm employs a direct solver and the other a pseudo-transient iterative solver. We will present results of our numerical simulations and discuss their application to natural shear zones.

Polymorphic transformations are key tracers of metamorphic processes, also used to estimate the pressure and temperature conditions reached by a rock. In particular, the quartz-coesite transition is commonly used to define the lower boundary of the ultrahigh-pressure (UHP) metamorphic field. The partial preservation of coesite included in garnets from UHP rocks bring considerable insights into the burial and exhumation mechanisms of the continental crust involved in convergent zone. Coesite was first described in the Western Alps by Chopin[1], in the Dora-Maria whiteschist, one of the most emblematic UHP rock worldwide. Although the partial preservation of coesite inclusions in garnet has long been attributed to the pressure vessel effect, the interrelationship and relative timing between fracturing and retrogression is still contentious.Here we study the reaction-deformation relationships of coesite inclusions initially enclosed in garnet and transforming into quartz during the decompression process. We combine 2D numerical thermo-mechanical models constrained by pressure-temperature-time (P-T-t) estimates from the Dora-Maira whiteschist. The model accounts for a compressible visco-elasto-plastic rheology including a pressure-density relationship of silica based on thermodynamic data. This allows us to study the effect of reaction-induced volume increase during decompression. Our results capture the typical fracture patterns of the host garnet radiating from retrogressed coesite inclusions and can be used to study the relative role of volume change associated with a change of P-T conditions on the style of deformation during decompression.The mechanisms of the coesite-quartz transformation and geodynamic implications are presented and validated against geological data. The effect of fluids on the phase transition and the conditions of access of fluids during the transformation are discussed in the light of the results of the thermo-mechanical models.This study demonstrates the high potential of thermo-mechanical modelling in enhancing our understanding of the processes involved in the formation and evolution of metamorphic minerals. [1]Chopin (1984) Contributions to Mineralogy and Petrology 86, 2, 107-118

Strain localisation is a key process that allows for the emergence of tectonic plates and controls their long-term deformation. Upper crustal levels are relatively cold and their rheology is thus governed by frictional plasticity. In order to predict the formation of tectonic plates and quantify the deformation of the Earth's upper shell, geodynamic modelling simulation tools must reliably account for deformation in the frictional plastic realm. Nevertheless, the simulation of frictional plastic strain localisation poses severe issues. Commonly employed implementations (visco-plastic and visco-elasto-plastic) often fail to accurately satisfy force balance and suffer from a lack of convergence upon mesh refinement. These problems are intimately linked to the fact that commonly employed models do not encompass any characteristic spatial or temporal scales of localisation. Various regularisation techniques can thus be used as a remedy. Here we investigate three popular regularisation techniques, namely viscoplasticity, gradient plasticity and the use of a Cosserat medium, and discuss their potential application for geodynamic modelling.

Julia(https://julialang.org) recently emerged as a very powerful high-level computer language for (parallel) scientific computing, which allows you to “write codes like in MATLAB”, while “achieving the speed of Fortran/C”. A particular strength of Julia is that it is easy to write composable software packages that talk to each other. Here we will discuss our efforts in making Julia a development platform for geodynamic applications that significantly simplifies the process of going from a working solver to a production code which runs on massively parallel (GPU) machines. We are working on a number of open-source packages that simplify certain steps that many geodynamics codes have in common:GeoParams.jl (https://github.com/JuliaGeodynamics/GeoParams.jl) is a package in which you can specify constitutive relationships (e.g., creeplaws). It automatically handles the (non-)dimensionalization of all input parameters, includes pre-defined creep laws (e.g., dislocation and diffusion creep laws), plotting routine and includes computational routines that can be directly integrated in your code. PETSc.jl (https://github.com/JuliaParallel/PETSc.jl) is the main interface from Julia to PETSc, including MPI support and automatic installations of PETSc (one of the main hurdles that existing users faced). We have recently extended the package to include an interface to DMSTAG, such that you create a staggered finite difference grid and assemble the stiffness matrix in a straightforward manner. You can use automatic differentiation tools in Julia to create the Jacobians for nonlinear equations, which again minimizes the required lines of code (compared to their C counterparts). At the same time, the full range of (nonlinear multigrid) PETSc solvers is available. This is thus very well suited to write implicit solvers. ParallelStencil.jl (https://github.com/omlins/ParallelStencil.jl) and ImplicitGlobalGrid.jl (https://github.com/eth-cscs/ImplicitGlobalGrid.jl) are packages that are devoted to solving stencils in a very efficient manner on (parallel) GPU or CPU machines, which scales to very large GPU-based computers. It is particularly efficient in combination with pseudo-transient iterative solvers and allow running codes on modern architectures. GeophysicalModelGenerator.jl (https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl) is a package that gives you a simple way to collect geophysical/geological data of a certain region and combine that to construct a 3D geodynamic input model setup. Ongoing efforts include the development of a grid generation and a marker and cell advection package that work, seamlessly with both ParallelStencil and PETSc. This will allow developers to apply both direct-iterative and pseudo-transient implicit solvers to the same problem, while only having to make minimal changes to the model setup. Combined, these packages will make the step from developing a new (nonlinear) solver to having an efficient (3D) production code much easier.

The development of highly efficient, robust and scalable numerical algorithms lags behind the rapid increase in massive parallelism of modern hardware. We address this challenge with the accelerated pseudo-transient iterative method and present here a physically motivated derivation. We analytically determine optimal iteration parameters for a variety of basic physical processes and confirm the validity of theoretical predictions with numerical experiments. We provide an efficient numerical implementation of pseudo-transient solvers on graphical processing units (GPUs) using the Julia language. We achieve a parallel efficiency over 96 % on 2197 GPUs in distributed memory parallelisation weak scaling benchmarks. 2197 GPUs allow for unprecedented terascale solutions of 3D variable viscosity Stokes flow on 49953 grid cells involving over 1.2 trillion degrees of freedom. We verify the robustness of the method by handling contrasts up to 9 orders of magnitude in material parameters such as viscosity, and arbitrary distribution of viscous inclusions for different flow configurations. Moreover, we show that this method is well suited to tackle strongly nonlinear problems such as shear-banding in a visco-elasto-plastic medium. A GPU-based implementation can outperform CPU-based direct-iterative solvers in terms of wall-time even at relatively low resolution. We additionally motivate the accessibility of the method by its conciseness, flexibility, physically motivated derivation and ease of implementation. This solution strategy has thus a great potential for future high-performance computing applications, and for paving the road to exascale in the geosciences and beyond., Geoscientific Model Development Discussions, ISSN:1991-962X, ISSN:1991-9611

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 W...

International audience; The preshortening Cretaceous Pyrenean Rift is an outstanding geological laboratory to investigate the effects of a pre‐rift salt layer at the sedimentary base on lithospheric rifting. The occurrence of a pre‐rift km‐scale layer of evaporites and shales promoted the activation of syn‐rift salt tectonics from the onset of rifting. The pre‐ and syn‐rift sediments are locally affected by high‐temperature metamorphism related to mantle ascent up to shallow depths during rifting. The thermo‐mechanical interaction between décollement along the pre‐existing salt layer and mantle ascent makes the Cretaceous Pyrenean Rifting drastically different from the type of rifting that shaped most Atlantic‐type passive margins where salt deposition is syn‐rift and gravity‐driven salt tectonics has been postrift. To unravel the dynamic evolution of the Cretaceous Pyrenean Rift, we carried out a set of numerical models of lithosphere‐scale extension, calibrated using the available geological constraints. Models are used to investigate the effects of a km‐scale pre‐rift salt layer, located at the sedimentary cover base, on the dynamics of rifting. Our results highlight the key role of the décollement layer at cover base that can alone explain both salt tectonics deformation style and high‐temperature metamorphism of the pre‐rift and syn‐rift sedimentary cover. On the other hand, in the absence of décollement, our model predicts symmetric necking of the lithosphere devoid of any structure and related thermal regime geologically relevant to the Pyrenean case.

International audience; The generation of ductile shear zones is essential for the formation of tectonic plate boundaries, such as subduction or strike-slip zones. However, the primary mechanism of ductile strain localization is still contentious. We study here the spontaneous generation of ductile shear zones by thermal softening using thermo-mechanical numerical simulations for linear and power-law viscous flow in one-dimension (1D), 2D and 3D. All models are velocity-driven. The 1D model exhibits bulk simple shear whereas the 2D and 3D models exhibit bulk pure shear. The initial conditions include a small temperature perturbation in otherwise homogeneous material. We use a series of 1D simulations to determine a new analytical formula which predicts the temperature evolution inside the shear zone. This temperature prediction requires knowledge of only the boundary velocity, flow law and thermal parameters, but no a priori information about the shear zone itself, such as thickness, stress and strain rate. The prediction is valid for 1D, 2D and 3D shear zones in bulk pure and simple shear. The results show that shear heating dominates over conductive cooling if the relative temperature increase is >50 °C. The temperature variation induced by the shear zone is nearly one order of magnitude wider than the corresponding finite strain variation so that no significant temperature variation occurs between shear zone and wall rock. Applying typical flow laws for lithospheric rocks shows that shear zone generation by thermal softening occurs for typical plate tectonic velocities of few cm.yr−1 or strain rates between and s−1. Shear stresses larger than 200 MPa can already cause strain localization. The results indicate that thermal softening is a feasible mechanism for spontaneous ductile shear zone generation in the lithosphere and may be one of the primary mechanisms of lithospheric strain localization.

Subduction zones are one of the main features of plate tectonics, they are essential for geochemical cycling and are often a key player during mountain building. However, several processes related to subduction zones remain elusive, such as the initiation of subduction or the exhumation of (ultra)high-pressure rocks.Collision orogens, such as the European Alps or Himalayas, provide valuable insight into long-term subduction zone processes and the larger geodynamic cycles of plate extension and subsequent convergence. For the Alps, geological reconstructions suggest a horizontally forced subduction initiation caused by the onset of convergence between the Adriatic and European plates. During Alpine orogeny, the Piemont-Liguria basin and the European passive magma-poor margin (including the Briançonnais domain) were subducted below Adria. The petrological rock record indicates burial and subsequent exhumation of both continental and oceanic crustal rocks that were exposed to (ultra)high-pressure metamorphic conditions during their Alpine burial-exhumation cycle. Moreover, estimates of exhumation velocities yield magnitudes in the range of several mm/yr to several cm/yr. However, published estimates of exhumation velocities, ages of peak metamorphic conditions and estimates for peak pressure and peak temperature vary partly significantly, even for the same tectonic unit. Consequently, many, partly significantly, contrasting conceptual models for subduction initiation (convergence versus buoyancy driven) or rock exhumation (channel-flow, diapirism, episodic regional extension, erosion etc.) have been proposed for the Alps. Complementary to the data-driven approach, mathematical models of the lithosphere and upper mantle system are useful tools to investigate geodynamic processes. These mathematical models integrate observational and experimental data with the fundamental laws of physics (e.g. conservation of mass, momentum and energy) and are useful to test conceptual models of subduction initiation and rock exhumation. Here, we present numerical solutions of two-dimensional petrological-thermo-mechanical models. The initial model configuration consists of an isostatically and thermally equilibrated lithosphere, which includes mechanical heterogeneities in the form of elliptical regions with different effective viscosity. We model a continuous geodynamic cycle of subsequent extension, no far-field deformation and convergence. During extension, the continental crust is necked, separated and mantle is exhumed, forming a marine basin bounded by passive margins. During the subsequent stage with no far-field deformation, the thermal field of the lithosphere is re-equilibrated above a convecting mantle. During convergence, subduction is initiated at one passive margin and the mantle lithosphere below the marine basin as well as the other passive margin are subducted. During progressive subduction, parts of the subducted continental upper crust are sheared-off the subducting plate and are exhumed to the surface, ultimately forming an orogenic wedge. For the convergence, we test the impact of serpentinites at the top of the exhumed mantle on orogenic wedge formation. We compare the model results with observational and experimental constraints, discuss the involved processes and driving forces and propose a model for subduction initiation and (ultra)high-pressure rock exhumation for the Alps.

Localization and softening mechanisms in a deforming lithosphere are important for subduction initiation or the generation of tectonic nappes during orogeny. Many localization mechanisms have been proposed as being important during the viscous, creeping, deformation of the lithosphere, such as thermal softening, grain size reduction, reaction-induced softening or anisotropy development. However, which localization mechanism is the controlling one and under which deformation conditions is still contentious. In this contribution, we focus on strain localization in viscous material due to the generation of anisotropy, for example due to the development of a foliation. We numerically model the generation and evolution of anisotropy during two-dimensional viscous deformation in order to quantify the impact of anisotropy development on strain localization and on the effective softening. We use a pseudo-transient finite difference (PTFD) method for the numerical solution. We calculate the finite strain ellipse during viscous deformation. The aspect ratio of the finite strain ellipse serves as proxy for the magnitude of anisotropy, which determines the ratio of normal to tangential viscosity. To track the orientation of the anisotropy during deformation, we apply the so-called director method. We will present first results of our numerical simulations and discuss their application to natural shear zones.

In the marker-in-cell method combined with staggered finite differences, Lagrangian markers carrying information on material properties are advected with the velocity field interpolated from the staggered Eulerian velocity grid. With existing schemes, velocity interpolation from the grid points to markers violates (to some extent) mass conservation requirement that causes excess convergence/divergence of markers and opening marker gaps after significant amount of advection. This effect is especially well visible in case of diagonal simple shear deformation along planes that are oriented at 45 degrees to the grid and marker circulation through grid corners.Here, we present a new second order velocity interpolation scheme that guarantees exact interpolation of normal strain rate components from pressure nodes (i.e. from the locations where these components are defined by solving of the mass conservation equation). This new interpolation scheme is thus applicable to both compressible and incompressible flow and is trivially expendable to 3D and to non-regular staggered grids.Performed tests show that, compared to other velocity interpolation approaches, the new scheme has superior performance in preserving continuity of the marker field during the long-term advection including the diagonal simple shear deformation and marker circulation through grid corners. We showcase a performance-oriented implementation of the new scheme using Julia language's shared memory parallelisation features. The Julia implementation of the new advection schemes further augments the ParallelStencil.jl related application collection with advection routines.

Ice’s predominantly viscous rheology exhibits a significant temperature and strain-rate dependence, commonly captured as a single deformation mechanism by Glen's flow law. However, Glen’s power-law relationship may fail to capture accurate stress levels at low and elevated strain-rates ultimately leading to velocity over- and under-estimates, respectively. Alternative more complex flow laws such as Goldsby rheology combine various creep mechanisms better accounting for micro-scale observations resulting in enhanced localisation of ice flow at glacier scales and internal sliding.The challenge in implementing Goldsby rheology arises with the need of computing an accurate partitioning of the total strain-rate among the active creep mechanisms. Some of these mechanisms exhibit grain-size evolution sensitivity potentially impacting the larger scale ice dynamics.We here present a consistent way to compute the effective viscosity of the ice using Goldsby rheology for temperature and strain-rate ranges relevant to ice flow. We implement a local iteration procedure to ensure accurate implicit partitioning of the total strain-rate among the active creep mechanisms including grain-size evolution. We discuss the composite deformation maps and compare the results against Glen's flow law. We incorporate our implicit rheology solver into an implicit 2D thermo-mechanical ice flow solver to investigate localisation of ice flow over variable topography and in shear margin configurations. We quantify discrepancies in surface velocity patterns when using Goldsby rheology instead of Glen's flow law.

We present two-dimensional thermo-mechanical numerical models of coupled lithosphere-mantle deformation, considering the upper mantle down to a depth of 660 km. We consider visco-elasto-plastic deformation and for the lithospheric and upper mantle a combination of diffusion, dislocation and Peierls creep. Mantle densities are calculated from petrological phase diagrams (Perple_X) for a Hawaiian pyrolite. The model generates a 120 Myrs long geodynamic cycle of subsequent extension (30 Myrs), cooling (70 Myrs) and convergence (20 Myrs) in a single and continuous simulation with explicitly modelling convection in the upper mantle. During lithosphere extension, the models generate an approximately 400 km wide basin of exhumed mantle bounded by hyper-extended passive margins. The models show that considering only the thermal effects of upper mantle convection by using an effective thermal conductivity generates results of lithosphere hyper-extension that are similar to the ones of models that explicitly model the convective flow. Applying a lower viscosity limit of 5 × 1020 Pa s suppresses convection and generates results different to the ones for simulations with a low viscosity asthenosphere having minimal viscosity of approximately 1019 Pa s. During cooling without far-field deformation, no subduction of the exhumed mantle is spontaneously initiated. Density differences between lithosphere and mantle are too small to generate a buoyancy force exceeding the mechanical strength of the lithosphere. The extension and cooling stages generate self-consistently a structural and thermal inheritance for the subsequent convergence stage. Convergence initiates subduction of the exhumed mantle at the transition to the hyper-extended margins. The main mechanism of subduction initiation is thermal softening for a plate driving force (per unit length) of approximately 15 TN m−1. If convection in the mantle is suppressed by high effective thermal conductivities or high, lower viscosity limits, then subduction initiates at both margins leading to divergent double-slab subduction. Convection in the mantle assists to generate a single-slab subduction at only one margin, likely due to mantle flow which exerts an additional suction force on the lithosphere. The first-order geodynamic processes simulated in the geodynamic cycle of subsequent extension, cooling and convergence are applicable to orogenies that resulted from the opening and closure of embryonic oceans bounded by magma-poor hyper-extended passive margins, which might have been the case for the Alpine orogeny.

The accommodation of motion on faults spans a large spectrum of slip modes, ranging from stable creep to earthquakes. While seismic slip modes certainly have the largest impact on the surface due to the induced ground shaking, it has been recognized that slow aseismic slip modes relax most of the accumulated stresses on a fault. It has also been suggested that aseismic slip controls seismic events, thus making this kind of slip mode key for earthquake prediction.Despite the importance of aseismic slow slip, its underlying physical mechanisms are still unclear. Commonly, slow slip events are modeled in terms of frictional failure, employing a rate-and-state model of fault friction, often also invoking fluids that alter frictional properties on the fault. However, at larger depths, frictional processes become increasingly difficult to activate due to the increase in ambient pressure and ductile processes are more likely to dominate deformation.Here we therefore investigate deep aseismic slip processes governed by ductile deformation mechanisms using 2D numerical models, where we employ a composite viscoelastic rheology combined with grain size reduction and shear heating as weakening processes. We show that the collaborative action of these two weakening mechanisms is sufficient to create the entire spectrum of aseismic slip, ranging from stable creep to long-term slow slip events. The results show that ductile deformation does not necessarily result in stable slip and induces slip modes with considerably larger velocities than the far-field plate velocities. Moreover, the propagation of ductile ruptures induces large stresses in front of the rupture tip which may also trigger short-term seismic events.

Although so-called "magma-poor" rifted margins display a large variability on a local scale, they are characterized by a number of common primary features worldwide such as their first-order architecture (proximal, necking, hyperextended, exhumation and oceanic domains), their lithological evolution along dip and the deformation processes associated with their different rifting stages. In this contribution, we first emphasize the primary morphological and lithological architecture of magma-poor rifted margins and how they relate to specific deformation modes (pure shear thinning, mechanical necking, frictional extensional wedge, detachment faulting and seafloor spreading). Second, we focus on the necking stage of rifting, which corresponds to the first major thinning event (when the crust is thinned from its initial thickness to ~ 10 km). We display the range of possible topographic and thermal evolutions of "magma-poor" and "sedimentary starved" rift systems depending on their lithosphere rheology. Our two-dimensional thermo-mechanical numerical models suggest that extension of lithospheres where the crust and the mantle are mechanically decoupled by a weak lower crust results in a complex morphotectonic evolution of rift systems, with formation of temporary restricted sub-basins framed by uplifted parts of the future distal margin. Mechanical decoupling between the crust and the mantle controls also largely the thermal evolution of rift systems during the necking phase since for equivalent extension rates and initial geotherms: (i) weak/decoupled lithospheres have a higher geothermal gradient at the end of the necking phase than strong/coupled lithospheres; and (ii) weak/decoupled lithospheres show intense heating of the lower crust at the rift center and intense cooling of the crust on either side of the rift center, unlike strong/coupled lithospheres. These behaviors contrast with the continuous subsidence and cooling predicted by the commonly used depth-uniform thinning model.

In this study, we use a state-of-the-art 2D numerical algorithm solving the standard thermo-mechanically coupled equations of continuum mechanics for slow flowing viscoelastoplastic material to model the evolution of rifting, thermal relaxation and convergence-to-collision of Alpine-type orogens in three stages. (1) A ca. 360 km wide basin that is floored by exhumed mantle and bounded by two conjugate magma-poor hyper-extended passive margins is generated during a 50 Myrs rifting period. An absolute extension velocity of 1 cm/yr is applied. (2) The passive margin system is thermally equilibrated during a subsequent cooling period of 60 Myrs without significant deformation in the lithosphere (no extension velocity). At this stage, we parameterise a serpentinization front on top of the exhumed mantle by replacing the dry peridotitic mantle by serpentinized mantle in one series of simulations. The thermally equilibrated system is used as a self-consistently generated initial configuration for the subsequent period of convergence lasting for 70 Myrs applying an absolute convergence velocity of 1.5 cm/yr. Values for the duration of deformation periods and for deformation velocities are chosen to allow for comparison between simulation results and petrological data from the Central and Western Alps. Density of all materials is either precomputed for characteristic bulk rock compositions and read in from precomputed thermodynamic look-up tables (Perple_X), or calculated during run time via a linearized equation of state (EOS). We quantify (1) the impact of a serpentinization front of the exhumed mantle on the subduction dynamics by increasing systematically the strength of the serpentinites, (2) the peak pressure and temperature conditions of subducted crustal material from the passive margins of the overriding and subducting plate by tracking pressure (P)-temperature (T)-time (t)-depth (z) paths of selected particles and (3) the driving forces of the system. Last, (4) the impact of metamorphic phase transitions is investigated by parameterising densification of crustal material. We compare the results of simulations in which density is computed as a simple linearized EOS to results of simulations in which density is a more realistic function of P and T using precomputed thermodynamic look-up tables. We discuss geometric similarities between the simulation results and 2D geodynamic reconstructions from field data, quantify the P-T-t-z-history of selected particles and compare it to P-T-t data obtained from natural rocks. First results indicate that the strength of the serpentinites controls whether the deformation within the orogenic core is driven by buoyancy forces (subduction channel model) or by far-field tectonic forces (orogenic wedge model). There is a transition from subduction channel to orogenic wedge model from low to intermediate strength of the serpentinites.

Shortening of Archaean and Paleoproterozoic continental lithospheres: large strains, but no orogeny Denis Gapais1, Jonathan Poh1, Philippe Yamato1, Thibault Duretz1, Florence Cagnard2 (1) Géosciences Rennes, UMR CNRS 6118, Université de Rennes 1, 35042 Rennes cedex, France (2) Bureau de Recherche Géologique et minière, 3 avenue Claude-Guillemin, BP 36009 45060 Orléans Cedex 2, France Denis.gapais@univ-rennes1.fr, jonathanpoh87@gmail.com, philippe.yamato@gmail.com, thibault.duretz@univ-rennes1.fr, f.cagnard@brgm.fr In many ancient deformation belts of Archaean and Paleoproterozoic age (e.g. Terre Adélie in East Antarctica, Finnish Svecofennides in Southern Finland, Murchison Belt in South Africa, Thompson Nickel Belt in Manitoba, Dharwar Craton in western India, Abitibi sub-Province in Québec, Trans-Hudson belt of Canada, Trans-Amazonian belt of Suriname), latest recorded deformations are compressive or transpressive. In these belts that involved hot and weak continental crusts, deformations are distributed with basically vertical tectonics and important crustal thickening. On the other hand, there is no evidence of syn-orogenic extension or late-orogenic collapse, as classically observed in modern orogens where extensional detachments are widespread.Analogue and numerical models emphasize that shortening of weak and hot lithospheres basically favour downward motions, which result in limited topographies. Field evidence further point to metamorphic isogrades rather parallel to the Earth surface at belt scale. Hence, metamorphic conditions are rather monotonous at the scale of individual belts, with limited metamorphic jumps and typical P-T paths with no significant adiabatic retrograde segments. Consistently, localized deep detrital sedimentary basins like foreland or intra-mountain basins, are not documented. Sedimentary records rather suggest distributed sedimentation processes. In addition, several lines of evidence tend to point out that cooling of ancient hot deformation belts was rather slow, which is consistent with distributed topographies and long-lasting erosion-driven exhumation processes. On these bases, we propose that gravity-driven collapse had no reason to occur in ancient hot deformation belts because important topographic gradients and orogeny could not develop as observed in modern mountain chains.

Numerous geological processes are governed by thermal and mechanical interactions. In particular, tectonic processes such as ductile strain localization can be induced by the intrinsic coupling that exists between deformation, energy and rheology. To investigate this thermomechanical feedback, we have designed 2-D codes that are based on an implicit finite-difference discretization. The direct-iterative method relies on a classical Newton iteration cycle and requires assembly of sparse matrices, while the pseudo-transient method uses pseudo-time integration and is matrix-free. We show that both methods are able to capture thermomechanical instabilities when applied to model thermally activated shear localization; they exhibit similar temporal evolution and deliver coherent results both in terms of nonlinear accuracy and conservativeness. The pseudo-transient method is an attractive alternative, since it can deliver similar accuracy to a standard direct-iterative method but is based on a much simpler algorithm and enables high-resolution simulations in 3-D. We systematically investigate the dimensionless parameters controlling 2-D shear localization and model shear zone propagation in 3-D using the pseudo-transient method. Code examples based on the pseudo-transient and direct-iterative methods are part of the M2Di routines (R¨ass et al., 2017) and can be downloaded from Bitbucket and the Swiss Geocomputing Centre website.

The model of an orogenic wedge has been applied to explain the tectonic evolution of many orogens worldwide. Orogenic wedges are characterized by (1) a first-order shear zone which underthrusts the mantle lithosphere and lower crust beneath the adjacent mantle lithosphere and (2) a sequence of second-order upper crustal shear zones which form tectonic nappes. Shear zone formation in the lithosphere is, however, incompletely understood. We perform two dimensional thermo-mechanical numerical simulations of lithospheric shortening to study shear zone formation, propagation and associated wedge formation. The only perturbation in the model lithosphere is a different temperature at the left (1300 ° C) and right (1400 ° C) half of the model bottom. Despite this smooth and weak perturbation, simulations show self-consistent and spontaneous formation of first- and second-order shear zones generating an orogenic wedge. The shear zones are caused by thermal softening and temperature-dependent rock viscosity. Lateral spacing of upper crustal shear zones spans between 30 and 50 km and is controlled by the depth of the boundary between upper and lower crust which acts as mechanical detachment level. Modelled upper crustal shear zones are active for ∼1 to ∼4 My. Surface processes such as sedimentation and erosion influence shear zone orientation, spacing and duration but do not impact fundamental processes of shear zone formation and propagation. Simulations produce both singly-vergent and doubly-vergent wedges. Topographic uplift rates are controlled by the applied bulk shortening rate. The modelled surface uplift and subsidence associated with crustal shear zones could explain major consecutive thrusting events and related sedimentation within flexural basins during the formation of the Helvetic nappe system in Western Switzerland. Furthermore, results for shear zone propagation along a mid-crustal detachment and associated uplift of crustal basement could potentially explain foreland basement-cored uplifts in natural orogens such as the Laramide orogen, Taiwan or the Shillong Plateau.

International audience; Studies of lithospheric mantle from (ultra)slow spreading ridges have shown that meltextraction at mid-ocean ridges may be incomplete, producing metasomatism/refertilizationof the shallow lithospheric mantle. However, it remains unclear whether similar processesoperate off axis and could affect the cooling lithosphere. Here, we report the discovery of agarnet xenocryst in a petit-spot lava sampled on the top of the downgoing Pacific plate infront of Japan. The trace-element composition of this garnet xenocryst, in particular thelow chromium, excludes a peridotitic origin, while the flat mid– to heavy rare earth elementpattern does not support direct crystallization from melt percolating through the oceaniclithosphere. Garnet formation is therefore interpreted as formed by a subsolidus reactionof a plagioclase-bearing cumulate during the progressive off-axis cooling of the lithosphere.Combining lithosphere cooling models and the specific physical conditions required for sub-solidus formation of garnet in tholeiitic systems (0.7–1.2 GPa) indicates that melt percolationto produce plagioclase-bearing cumulate occurs >150 km off axis. These conditions supportthat some low-degree melts produced off axis are not collecting to form mid-oceanic ridgebasalt (MORB), but percolate and crystallize during the cooling and thickening of the litho-spheric mantle. The demonstration of mantle refertilization affecting Pacific lithosphere offaxis is critical because such a process could explain the presence of metasomatic domainswith distinct physical and chemical properties in the depleted oceanic lithosphere.

Natural pinch-and-swell structure in centimetre-thick calcite layers shows a reduction of grain size from swell towards pinch. However, the impact of grain size evolution on necking and pinch-and-swell formation is incompletely understood. We perform zero-dimensional (0D) and 2D thermo-mechanical numerical simulations of calcite layer extension to quantify the impact of grain size evolution on necking for bulk extension rates between 10 −12 s −1 and 10 −14 s −1 and temperatures around 350 °C. For a combination of diffusion and dislocation creep we calculate grain size evolution according to the paleowattmeter (grain size is proportional to mechanical work rate) or the paleopiezometer (grain size is proportional to stress). Numerical results fit three observations: (i) significant thickness variations along the layer after extension, (ii) grain size reduction from swells towards pinches, and (ii) dislocation creep dominated deformation in swells and significant contribution of diffusion creep in pinches. Modelled grain size in pinches (10–60 μm) and swells (70–800 μm) is close to observed grain size in pinches (21 ± 6 μm) and in swells (250–1500 μm). In the models, grain size evolution has a minor impact on necking, and viscous shear heating and grain size evolution have a negligible thermal impact.

Recent development of many multiphysics modeling tools reflects the currently growing interest for studying coupled processes in Earth Sciences. The core of such tools should rely on fast and robust mechanical solvers. Here we provide M2Di, a set of routines for 2-D linear and power law incompressible viscous flow based on Finite Difference discretizations. The 2-D codes are written in a concise vectorized MATLAB fashion and can achieve a time to solution of 22 s for linear viscous flow on 10002 grid points using a standard personal computer. We provide application examples spanning from finely resolved crystal-melt dynamics, deformation of heterogeneous power law viscous fluids to instantaneous models of mantle flow in cylindrical coordinates. The routines are validated against analytical solution for linear viscous flow with highly variable viscosity and compared against analytical and numerical solutions of power law viscous folding and necking. In the power law case, both Picard and Newton iterations schemes are implemented. For linear Stokes flow and Picard linearization, the discretization results in symmetric positive-definite matrix operators on Cartesian grids with either regular or variable grid spacing allowing for an optimized solving procedure. For Newton linearization, the matrix operator is no longer symmetric and an adequate solving procedure is provided. The reported performance of linear and power law Stokes flow is finally analyzed in terms of wall time. All MATLAB codes are provided and can readily be used for educational as well as research purposes. The M2Di routines are available from Bitbucket and the University of Lausanne Scientific Computing Group website, and are also supplementary material to this article.

Tectonic nappes are observed for more than a hundred years. Although geological studies often refer to a nappe theory, the physical mechanisms of nappe formation are still incompletely understood. We apply two-dimensional numerical simulations of shortening of a passive margin, to investigate the thermo-mechanical processes of detachment, transport and stacking of nappes. We use a visco-elasto-plastic model with standard creep flow laws and Drucker-Prager yield criterion. We consider tectonic inheritance with two initial mechanical heterogeneities: (1) lateral heterogeneity of the basement-cover interface due to half-grabens and horsts and (2) vertical heterogeneities due to layering of mechanically strong and weak sedimentary units. The model shows detachment and horizontal transport of a thrust nappe and stacking of this thrust nappe above a fold nappe. The detachment of the thrust sheet is triggered by stress concentrations around the sediment-basement contact and the resulting brittle-plastic shear band formation. The horizontal transport is facilitated by a basal shear zone just above the basement-cover contact, composed of thin, weak sediments. Fold nappe formation occurs by a dominantly ductile closure of a half-graben and the associated extrusion of the half-graben fill. We apply our model to the Helvetic nappe system in Western Switzerland, which is characterized by stacking of the Wildhorn thrust nappe above the Morcles fold nappe. The modeled structures and temperature field agree with data from the Helvetic nappe system. The mechanical heterogeneities must generate contrasts in effective viscosity (i.e. ratio of stress to strain rate) of four orders of magnitude to model nappe structures similar to the ones of the Helvetic nappe system.

International audience; It is widely accepted that subduction initiation at modern Earth passive margin systems critically depends on the buoyancy and strength of the oceanic lithosphere and requires failure of the load-bearing crustal and mantle layers. As such, subduction initiation upon orthogonal convergence is controlled by the age of the oceanic lithosphere and thus the strength contrast at the margin. However, it is still unclear where along the margin and how subduction initiates. In particular, rheologically-controlled mechanisms are poorly understood and require further investigation. Therefore, this combined analogue and numerical modelling study aims at exploring the effects of first order rheological and kinematic conditions on subduction initiation at passive margins. Our results highlight the sensitivity of early stages of subduction initiation to the initial rheological setup as well as evolving thermo-mechanical feedback mechanisms. Additionally, they provide more insights on the conditions responsible for the locus of subduction initiation. We infer that the locus of subduction is controlled by the rheology of both the crust and the mantle lithosphere at the margin, which is in strong correlation with the thermal age of the oceanic lithosphere. Suitable conditions for subduction initiation at passive margins correspond to an intermediate age (ca. 80 Myr) of the oceanic lithosphere. In all other cases deformation localizes within the oceanic lithosphere or affects the entire continent. We advocate the significance of crust-mantle decoupling at the passive margin for its inversion and possible evolution towards a subduction zone. Additionally, the development of a self-sustaining subduction zone will fail in the absence of weakening mechanisms taking place in the mantle lithosphere such as thermal softening (shear heating) and low temperature plasticity (Peierls mechanism).

International audience; This article points out for the first time a striking correlation between the paleogeography of Upper Triassic deposits and the mode of crustal stretching around and inside the Northern Iberia plate during the Cretaceous transtensional event. First, we propose a review of the architecture of the basins, which opened during the mid-Cretaceous times along the Iberia-Eurasia plate boundary. Similarly to the emblematic Parentis basin, these basins exhibit a peculiar synclinal-shaped profile and are devoid of prominent block faulting. The top of the basement is characterized by gentle slopes, which dip symmetrically toward the center of the basins. As revealed by recent comparisons with geologically-constrained rifting models established from the North Pyrenean Zone, this architecture results from the thinning of a heterogeneous continental crust under greenschist facies conditions. Basement deformations are thus dominantly ductile and are characterized by large-scale boudinage, hyper-thinning and subsequent lateral extraction. Bulk deformation of the upper and middle crust leads to the formation of anastomosed shear zones and the development of mylonitic fabric. Tectonic lenses consisting of crustal material remain welded on top of the exhuming mantle. The common character shared by all the pre-rift sequences of the studied basins is the presence of a thick low-strength Upper Triassic evaporites and clays layer belonging to the Keuper group (i.e. pre-rift salt and clay unit). In the studied basins, efficient décollement along the pre-rift salt and clay unit triggers mechanical decoupling and gliding of the pre-rift cover that remains in the center of the basin. Thus, during the early rifting phase, the basement undergoes thinning while the pre-rift cover remains preserved in the basin center. In response to hyper-thinning and horizontal extraction of the continental crust, hot mantle is exhumed beneath the pre- and syn-rift cover. Subsequent thermal exchange (i) promotes ductile deformation of the basement and (ii) induces the development of HT-LP metamorphic conditions in the pre-rift sediments and at the base of the syn-rift flysch levels. This thermal event is well recorded in the axial portion of the Pyrenean realm (future North Pyrenean Zone) as well as in the pre-rift sediments of the Cameros basin (northern Spain). The architecture of the smooth-slopes type basins thus contrasts with the structure of Iberia-Newfoundland Atlantic margins which are characterized by (i) top-basement detachment faults accommodating crustal extension through rotation and translation of undeformed basement blocks, and (ii) the individualization of continental extensional allochthons tectonically emplaced over exhumed lower crust or mantle rocks. Finally, using recent paleogeographic reconstructions, we show that the distribution of the pre-rift salt and clay unit remarkably matches the distribution of the Pyrenean and peri-Pyrenean smooth-slopes type basins. This allows for the first time to propose a genetic link between the distribution of evaporite-bearing pre-rift sedimentary formations and the development of smooth-slopes basins.

We present two-dimensional numerical simulations for shortening a viscoelastoplastic lithosphere to quantify the impact of elasticity on strain localization due to thermal softening. The model conserves energy and mechanical work is converted into heat or stored as elastic strain energy. For a shear modulus G = 1010 Pa, a prominent lithospheric shear zone forms and elastic energy release increases the localization intensity (strain rate amplification). For G = 5 × 1010 Pa shear zones still form but deformation is less localized. For G = 1012 Pa, the lithosphere behaves effectively viscoplastic and no shear zones form during homogeneous thickening. Maximal shearing-related increase of surface heat flux is 15–25 mW m−2 and of temperature at lower crustal depth is ∼150 °C, whereby these peak values are transient (0.1–1 My). Elasticity and related energy release can significantly contribute to strain localization and plate-like behaviour of the lithosphere required for plate tectonics.

International audience; We present two-dimensional numerical models of plate subduction and collision inspired by the Scandinavian Caledonian orogeny to investigate the possible impact of continental crust partial melting on the exhumation of ultra-high pressure metamorphic rocks. Three possible reactions were tested: low temperature solidus representing H2O-fluid-saturated partial melting, and two end-member reaction curves for dehydration melting. Thermo-mechanical effects of partial melting were implemented as (1) a viscosity decrease as a determined rheologically critical melt percentage was reached (here 0.1), (2) a change in effective heat capacity and adiabatic heating/cooling accounting for a latent heat term in the heat equation. Among the 3 tested reactions, only H2O-fluid-saturated partial melting drastically modifies the collision dynamics from the non-melting reference model holding all other parameters constant. A substantially low general viscosity truncation (here 1017 Pa s) is needed to properly resolve the effect of partial melting on deep collision processes. Low temperature melting indeed induces the development of a low viscosity buoyant plume prior to slab detachment, where migmatites exhume from UHP conditions at rates and with pressure–temperature paths similar to the natural values acknowledged for the Norwegian Caledonides. High temperature melting has no drastic influence on early collision dynamics. While positive buoyancy remains the first order driver for the exhumation of buried continental rocks, exhumation initiates in these cases with eduction subsequent to slab detachment. Melting and formation of a migmatite plume can later occur along decompression path while continental crust undergoes thermal reequilibration at temperatures above 900 °C. Some of the partially molten material can also relaminate in the overriding plate rather than exhume within the collision zone. Even if minor in terms of amount of magma produced, H2O-fluid-saturated partial melting at UHP conditions could therefore have a dramatic rheological effect and actually limits continental rocks subduction and facilitates their exhumation.

International audience; The dynamics of magma flow is highly affected by the presence of a crystalline load. During magma ascent, it has been demonstrated that crystal-melt segregation constitutes a viable mechanism for magmatic differentiation. Moreover, crystal-melt segregation during magma transport has important implications in term of magma rheology, but also in term of differentiation of the continental crust. However, the influences of the crystal volume percentage (φ), of their geometry, their size and their density on crystal-melt segregation are still not well constrained. To address these issues, we performed a parametric study using 2D direct numerical simulations, which model the ascension of a crystal-bearing magma in a vertical dyke. Using these models, we have characterised the amount of segregation as a function of different physical properties including φ, the density contrast between crystals and the melt phase (Δρ), the size of the crystals (Ac) and their aspect ratio (R). Results show that small values of R do not affect the segregation. In this case, the amount of segregation depends upon four parameters. Segregation is highest when Δρ and Ac are large, and lowest for large pressure gradient (Pd) and/or large values of dyke width (Wd). These four parameters can be combined into a single one, the Snumber, which can be used to quantify the amount of segregation occurring during magma ascent. Based on systematic numerical modelling and dimensional analysis, we provide a first order scaling law which allows quantification of the segregation for an arbitrary Snumber and φ, encompassing a wide range of typical parameters encountered in terrestrial magmatic systems. Although developed in a simplified system, this study has strong implications regarding our understanding of crystal segregation processes during magma transport. Our first order scaling law allows to immediately determine the amount of crystal-melt segregation occurring in any given magmatic dyke system.

We investigate the dynamics of thermally activated shear localization in power law viscoelastic materials. A two-dimensional (2D) thermomechanical numerical model is applied that uses experimentally derived flow laws for rock. We consider viscous and viscoelastic rheologies and show that the numerical solutions for shear bands are mesh insensitive and energetically conservative. Deformation under long-term tectonic strain rates () yields to shear localization on the scale of kilometres. Although viscous and viscoelastic models exhibit comparable shear zone thickness, the timing of shear localization and stress magnitudes are affected by viscoelasticity. Large values of shear modulus ( Pa) promotes fast stress loading (within strain) and localization ( strain), whereas lower values ( Pa) delays localization ( strain) and reduces the maximum effective stress by a factor two. High stress exponents (up to ) focuses the deformation into narrow shear zones ( m) while maintaining a relatively high average stress...

We present the results of two-dimensional thermo-mechanical numerical simulations of the formation of kilometre-scale shear zones and tectonic nappes which are generated by thermal softening during horizontal shortening. This model is thermo-mechanically coupled and conserves energy by converting mechanical work into heat. The results show that in homogeneous material, the temperature increase in shear zones caused by viscous heating can be ~100 °C. The characteristic width of the thermally perturbed part of the deforming zone is three times larger than the width of the corresponding shear zone because of diffusive heating of the wall rock. Therefore, temperature variations initiated by shear heating may be difficult to identify based on petrological inferences of field data. During lithospheric shortening with a bulk rate of 2.5 × 10−15 s−1, a nappe forms in the crust at a depth of ~45 km. Peak temperatures in the nappe reach ~600 °C (~100 °C increase due to shear heating). Significant shear heating occurs only locally and is not active continuously in time along the entire crustal shear zone. Peak values of the maximal principal stress locally reach 2.5 GPa in the nappe. Peak values of temperature and maximal principal stress increase locally to ~770 °C (~300 °C increase due to shear heating) and ~3.8 GPa, respectively, if the applied bulk shortening rate is 5 × 10−15 s−1 and if higher effective viscosities for the crust are applied. In the model lithosphere, high differential stresses (>500 MPa) and significant stress deviations from the lithostatic pressure (>50%) occur only locally and are transient. The highest average effective viscosity of the model lithosphere is ~2 × 1022 Pa s, which is in agreement with the estimates for the continental lithosphere. A characteristic feature of the presented dynamic nappe model is a significant stress decrease and cooling within a hundred thousand to 2 million years. These rates of decompression and cooling are discussed in the context of microstructural and metamorphic observations. The results also show that differential stresses in mature shear zones can be several hundred MPa smaller than at the onset of shear zone formation, which suggests that estimates based on microstructural analysis of rocks in ductile shear zones provide minimum values. The presented model could explain the formation of coherent Alpine nappes with inferred high-pressure (1.5–2.5 GPa) and ultrahigh-pressure (>3 GPa) conditions in 40–60 km depth at the base of an orogenic wedge.

We study the formation of localized shear zones during the layer-parallel extension of viscous multi-layers using two-dimensional numerical simulations based on the finite-difference method. For power-law viscous layers and a linear viscous embedding medium, the extended multi-layer develops boudins due to necking. For power-law viscous layers embedded in a power-law viscous medium, the extended multi-layer develops first distributed necks, and subsequently a localized shear zone with a vertical offset (with a size of several layer thicknesses) along the multi-layer. During the extension, the deformation style switches from distributed and symmetric necking to localized and asymmetric shearing. A localized shear zone develops in the viscous multi-layer although the rheology is everywhere strain-rate-hardening (power-law stress exponent >1) and no material softening and/or energy feedback mechanism (e.g., shear heating combined with a temperature-dependent viscosity) is applied. The shear localization is caused by structural softening because the formation of a localized shear zone decreases the bulk resistance and hence the work required to deform the multi-layer. A localized shear zone forms in the multi-layer when the distance between the stiff layers is approximately equal to or less than the layer thickness. The shear localization was observed in multi-layers with nine and with only three stiff layers.

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.2 GPa with corresponding deviatoric stresses ~ 1.5 GPa and occur in a depth of ~ 42 km. Maximal pressure of ~ 3.4 GPa 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

We evaluate the parameters that control the thickness of ductile shear zones that are generated by viscous heating. We employ two-dimensional thermomechanical numerical models to simulate shear zone development under compression. Results show that the shear zone thickness is essentially independent on the numerical resolution and the initial size of a weak inclusion that triggers shear localization. A simple scaling law is derived which predicts the thickness with six physical parameters: far-field stress and strain rate, thermal conductivity (both constant and temperature dependent), initial temperature, activation energy, and stress exponent. The scaling law is confirmed by numerical simulations for a wide range of parameters. For crustal deformation conditions typical temperature increase ranges between 50 ◦ C and 200 ◦ C, and the predicted thickness is between 1 km and 7 km. These thicknesses agree with natural crustal and lithospheric shear zones suggesting that shear heating is a process controlling crustal shear zone formation.

We study the three-dimensional (3-D) deformation during detachment of a lithospheric slab with simple numerical models using the finite element method. An initially vertical layer of power law viscous fluid mimics the slab and is surrounded by a linear or power law viscous fluid representing asthenospheric mantle. We quantify the impact of slab size and shape (symmetric/asymmetric) on slab detachment and identify two processes that control the lateral (i.e., along-trench) slab deformation: (1) the horizontal deflection of the lateral, vertical slab sides (> 100 km with velocities up to 16 mm/yr) and (2) the propagation of localized thinning (necking) inside the slab (with velocities >9 cm/yr). The lateral propagation velocity is approximately constant during slab detachment. Larger slabs (here wider than approximately 300 km) detach with rates similar to those predicted by 2-D models, whereas smaller slabs detach slower. Implications for geodynamic processes and interpretations of seismic tomography are discussed.

The eastern Tibetan margin is characterized by a steep topographic gradient and remarkably lateral variations in crustal/lithospheric structure and thermal state. GPS measurements show that the surface convergence rate in this area is strikingly low. How can such a mountain range grow without significant upper crustal shortening? In order to investigate the formation mechanism of the eastern Tibetan-type margins, we conducted 2D numerical simulations based on finite difference and marker-in-cell techniques. The numerical models were constrained with geological and geophysical observations in the eastern Tibetan margin. Several major parameters responsible for topography building, such as the convergence rate, the erosion/sediment rate, and the presence of partially molten crust, were systematically examined. The results indicate that the presence of partially molten material in the middle/lower crust can make a positive contribution to the formation of steep topography, but it is not a necessary factor. A steep topographic gradient may be a characteristic feature when a thin lithosphere with thick crust converges with a thick lithosphere with thin crust. In the context of a high erosion rate, the Longmen Shan range still gains and maintains its steep high topography to the present. This could be explained by exerting a large push force on Tibet side. Our numerical experiments suggest that topographic characteristic across the eastern Tibetan-type margins is mainly derived from isostatic equilibration forces and intensive convergence between two continental lithospheres with totally different rheological properties.

[1] Here we describe a new staggered grid formulation for discretizing incompressible Stokes flow which has been specifically designed for use on adaptive quadtree-type meshes. The key to our new adaptive staggered grid (ASG) stencil is in the form of the stress-conservative finite difference constraints which are enforced at the “hanging” velocity nodes between resolution transitions within the mesh. The new ASG discretization maintains a compact stencil, thus preserving the sparsity within the matrix which both minimizes the computational cost and enables the discrete system to be efficiently solved via sparse direct factorizations or iterative methods. We demonstrate numerically that the ASG stencil (1) is stable and does not produce spurious pressure oscillations across regions of grid refinement, which intersect discontinuous viscosity structures, and (2) possesses the same order of accuracy as the classical nonadaptive staggered grid discretization. Several pragmatic error indicators that are used to drive adaptivity are introduced in order to demonstrate the superior performance of the ASG stencil over traditional nonadaptive grid approaches. Furthermore, to demonstrate the potential of this new methodology, we present geodynamic examples of both lithospheric and planetary scales models.

Common features visible on a large majority of the seismic cross sections are chimneys or pipe structures. They represent regions of focused fluid flow in porous media. As seismic surveys are widely performed in many regions where the subsurface is of economic interest, a better understanding of the formation and evolution process of these chimneys is vital. They should be considered when performing risk assessment linked to leakage within subsurface waste storage projects. They might also lead to a better understanding of fluid migration pathways and subsurface fluid localization. In that context, we propose a new physical model that predicts the formation and the evolution in space and time of these chimneys. We use a two-dimensional (2D) implicit solver and the three-dimensional (3D) high-resolution iterative parallel GPU code to solve a thermodynamically consistent system of nonlinear equations for two-phase flow in deforming porous media. We will show that the different 2D implicit and 3D iterative methods used to solve the fully coupled system of nonlinear equations are in good agreement. They predict the formation and the propagation of a nonlinear solitary wave, resulting in a chimney formation. These chimneys, with order of magnitude permeability increase, are a natural outcome of an interplay between buoyancy forces driving upward fluid propagation and a resistance of a deformable rock to locally increasing fluid pressure. We will also discuss and highlight the importance of a proper coupling between the geomechanics (Stokes solver) and the reservoir fluid flow (nonlinear Darcy solver).

SUMMARY Calculating surface topography in geodynamic models is a common numerical problem. Besides other approaches, the so-called ‘sticky air’ approach has gained interest as a free-surface proxy at the top boundary. The often used free slip condition is thereby vertically extended by introducing a low density, low viscosityfluid layer. This allows the air/crust interface to behave in a similar manner to a true free surface. We present here a theoretical analysis that provides the physical conditions under which the sticky air approach is a valid approximation of a true free surface. Two cases are evaluated that characterize the evolution of topography on different timescales: (1) isostatic relaxation of a cosine perturbation and (2) topography changes above a rising plume. We quantitatively compare topographies calculated by six different numerical codes (using finite difference and finite element techniques) using three different topography calculation methods: (i) direct calculation of topography from normal stress, (ii) body-fitting methods allowing for meshing the topography and (iii) Lagrangian tracking of the topography on an Eulerian grid. It is found that the sticky air approach works well as long as the term (ηst/ηch)/(hst/L) 3 is sufficiently small, where ηst and hst are the viscosity and thickness of the sticky air layer, and ηch and L are the characteristic viscosity and length scale of the model, respectively. Spurious lateral fluctuations of topography, as observed in some marker-based sticky air approaches, may effectively be damped by an anisotropic distribution of markers with a higher number of markers per element in the vertical than in the horizontal direction.