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

International audience; The nucleation and rupture processes of deep-focus earthquakes have remained enigmatic ever since their discovery. These earthquakes occur mostly within the mantle transition zone where brittle failure is extremely unlikely due to the elevated pressures at these depths. Hence, other mechanisms have to be invoked to explain the occurrence of these events. To date, two main hypotheses have been put forward to explain deep focus earthquakes: transformational faulting (due to the polymorphic phase change of metastable olivine to either wadsleyite or ringwoodite) and thermal runaway (due to the conversion of deformational work to heat). More recently, it has been proposed that the feedback between those two mechanisms may explain the observed two-stage ruptures of large deep-focus earthquakes.To better understand the potential feedback between transformational faulting and thermal runaway, it is necessary to determine the stresses induced by the phase change due to i) the grain size reduction and corresponding viscosity reduction of the transformed material and ii) the volume reduction of the transformed phase. The former process triggers a stress transfer from the transformed material to the untransformed material, whereas the latter results in elevated stresses around the transformed phase.In this study, we employ numerical models with a viscoelastic compressible rheology to quantify the stress levels and patterns resulting from both processes. To gain a better understanding of the parameters controlling the stress transfer from transforming regions to the surrounding matrix, we employ simplified numerical models where transforming regions are approximated using elliptical inclusions. In a second step, more realistic model geometries are used to additionally study the effect of the morphology of transformed regions on stress levels and heterogeneities.Results show that both processes result in significantly different stress evolution upon a phase transition. Whereas a phase transition affecting only the viscosity of the transformed material results in moderate stress increases which occur on relatively long timescales, a phase transition affecting both viscosity and density results in significantly larger stresses, which also exhibit a significantly faster build-up. In both cases, the attained stress are sufficiently large to activate additional ductile weakening mechanisms that could trigger ductile ruptures. The higher stress levels resulting from the combined effect of a viscosity and density change would likely result in stronger weakening effects and faster occurrence of ductile failure

Quantitative phase petrology analysis of two samples showing variable strain record and different peak pressure (P) conditions show that H2O is a key parameter in the development of contrasting high-pressure (HP) peak conditions in adjacent polymetamorphic rocks. A shear zone in a paragneiss from the Gran Paradiso nappe, Western Alps, shows an Alpine foliation and records a peak P of 1.9 GPa, for a temperature of 500 to 520 °C. A few tens of meters away from the shear zone, a paragneiss showing no apparent Alpine age deformation records a peak P of maximum 1.4 GPa for the temperature range of 500 to 520 °C. The H2O content of the latter has potentially been reduced to low contents following the pre-Alpine, Variscan amphibolite facies, and the absence of re-hydration prior to Alpine orogeny could have inhibited the formation of HP mineral assemblages. The validity of this interpretation is questioned here by considering the mechanical effect of H2O undersaturated rocks on deformation and P during deformation. Based on a thermo-mechanical numerical modelling study, we show that heterogeneities in fluid saturation conditions between rocks lead to strength contrasts that are sufficient to trigger a dynamic P in the range of several hundreds of MPa. In particular, the models successfully reproduce the measured peak P between the two paragneiss studied. This model could be applied to other H2O deficient rocks from HP tectonic units to further explore the role of H2O on the rheology and hence assess its potential impact in the preservation of low P bodies in otherwise HP units from continental collision settings.

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

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

International audience; Deformations of the colder regions of the lithosphere mainly occur in the frictional regime. In geodynamic models, frictional plastic deformations are often highly localised (shear bands) and are used as proxies for faults. However, capturing the generation and evolution of shear bands in geodynamic models is troublesome. Indeed, mesh dependency and lack of convergence affect, to some extent, the results of geodynamic models. Here we extend the most common plasticity implementation used in geodynamic codes (effective viscosity approach) to include the combined effects of elasto-plastic compressibility, plastic dilatancy, strain softening and viscoplasticity. The latter acts as a regularisation that cures most of the known issues of geodynamic models related to frictional plasticity. Using regularised models based on the M2Di MATLAB routines, we show that volumetric elasto-plastic deformations can significantly impact crustal-scale shear banding. We also show that the artificial overstress caused by viscoplasticity can be mitigated by employing power-law models. Furthermore, we demonstrate that plasticity algorithms common in geodynamics (based on the effective viscosity approach) can be as accurate as those obtained with algorithms typically used in engineering (return mapping with a consistent tangent operator). Finally, we show examples of long-term tectonic deformations using the state-of-the art geodynamic code MDoodz. They indicate that viscoplastic regularisation can be used efficiently to obtain reliable simulations in geodynamics.

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.

Ultrahigh-pressure (UHP) continental crustal rocks were first discovered in the Western Alps in 1984 and have since then been observed at many convergent plate boundaries worldwide. Unveiling the processes leading to the formation and exhumation of (U)HP metamorphic crustal rocks is key to understand the geodynamic evolution of orogens such as the Alps. Previous numerical studies investigating (U)HP rock exhumation in the Alps predicted deep (>80 km) subduction of crustal rocks and rapid buoyancy-driven exhumation of mainly incoherent (U)HP units, involving significant tectonic mixing forming so-called mélanges. Furthermore, these predictions often rely on excessive erosion or periods of divergent plate motion as important exhumation mechanism. Inconsistent with field observations and natural data, application of these models to the Western Alps was recently criticised. Here, we present models with continuous plate convergence, which exhibit local tectonic-driven upper plate extension enabling compressive- and buoyancy-driven exhumation of coherent (U)HP units along the subduction interface, involving feasible erosion. The two-dimensional petrological-thermo-mechanical numerical models presented here predict both subduction initiation and serpentinite channel formation without any a priori prescription of these two features. The (U)HP units are exhumed coherently, without significant internal deformation. Modelled pressure and temperature trajectories and exhumation velocities of selected crustal units agree with estimates for the Western Alps. The presented models support previous hypotheses of synconvergent exhumation, but do not rely on excessive erosion or divergent plate motion. Thus, our predictions provide new insights into processes leading to the exhumation of coherent (U)HP crustal units consistent with observations and natural data from the Western Alps.

International audience; We document the role of sedimentary burial and salt tectonics in controlling the deformation style of continental crust during hyperextension. The Iberian-European boundary records a complex history of Cretaceous continental extension, which has led to the development of so-called smooth-slope type basins. Based on the review of the available geological constraints (crustal-balanced cross sections, sedimentary profile evolution, RSCM thermometer, low-temperature thermochronology) and geophysical data (Bouguer anomaly, Moho depth, seismic reflection profiles, and Vp/Vs velocity models) on the Tartas, Arzacq, Cameros, Parentis, Columbrets, Mauléon, Basque-Cantabrian and Internal Metamorphic Zone basins, we shed light on the main characteristics of this type of basin. This synthesis indicates that crustal thinning was influenced by two decoupling horizons: the middle crust and Triassic pre-rift salt, initially located between the basement and pre-rift sedimentary cover. These two horizons remained active throughout basin formation and were responsible for depth-dependent thinning of the crust and syn-rift salt tectonics. We therefore identify several successive deformation phases involving (1) pure shear dominated thinning, (2) simple shear dominated thinning and (3) continental breakup. In the first phase, distributed deformation resulted in the development of a symmetric basin. Field observations indicate that the middle and lower crust were under dominantly ductile conditions at this stage. In the second phase, deformation was localized along a crustal detachment rooted between the crust and the mantle and connecting upwards with Triassic pre-rift salt. During continental breakup, basin shoulders recorded the occurrence of brittle deformation while the hyperextended domain remained under predominantly ductile thinning. The formation of smooth-slope type extensional basins was intrinsically linked to the combined deposition of thick syn-rift and breakup sequences, and regional salt tectonics. They induced significant burial and allowed the continental crust and the pre-rift sequence to deform under high temperature conditions from the rifting to continental breakup stages.

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.

International audience; Orogens are traditionally classified according to their tectonic style. Paleoproterozoic tectonics is referred to as “ancient-style tectonics” while Proterozoic tectonics is referred to as “modern-style tectonics”. Ancient-style tectonics is characterised by distributed vertical structures and low topography gradients, often associated with diapirism and partial melting. In contrast, modern-style tectonics involve prominent strain localisation and the formation of thrusts, nappes and high topographic gradients. However, the parameters controlling the transition from ancient to modern-style tectonics are poorly understood. To quantify this transition, a combination of 1D and 2D high resolution lithospheric-scale thermo-mechanical models was conducted. The parameters controlling the strength of the lithosphere (i.e., Moho temperature, strain rate, crustal rheology, crustal radiogenic heat production and role of shear heating) were investigated in detail. Our results show that tectonic style is controlled by the maximum of crustal strength (shear stress). Modern-style tectonics is observed to occur when the maximum of crustal strength is greater than 300 MPa. At the opposite, a maximum crustal strength lower than 300 MPa leads to ancient style tectonic structures. Therefore, crustal rheology, temperature and background strain rate significantly influence the transition from ancient to modern-style tectonics. Shear heating remains a key factor in promoting strain localisation in modern-style tectonics. Crustal radiogenic heat production has a moderate influence by increasing/decreasing the tendency for faulting within the crust. This crustal strength criterion also provides an excellent fit for a second potential proxy: a localisation criterion of ca. 225°C. These two proposed proxies can be used interchangeably to predict the transition from ancient to modern-style tectonics.

International audience; We used high-resolution (500 × 250 m) two-dimensional lithospheric-scale thermo-mechanical numerical modeling to unravel the unexpected topographic and thermal evolutions recorded during the necking phase of several rift systems worldwide. Through a systematic analysis we studied how the lithosphere rheology impacts the topographic and thermal evolutions across the entire width of magma-poor and sediment-starved rift systems until their crust is locally thinned to 10 km. We quantified the evolution of topography, uplift and subsidence rates, accommodation and emerged space creation, temperature and surface heat flow for a wide panel of crustal and mantle rheologies to provide an overview of possible rifting behaviors. Extension of a lithosphere for which the crust and mantle are mechanically decoupled by a weak lower crust generates complex morphotectonic evolutions, with the formation of temporarily restricted sub-basins framed by uplifted parts of the future distal margin. Mechanical decoupling between the crust and 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. Accommodation space in the evolving basins is first generated by vertical crustal velocities and subsequently by horizontal velocities causing the widening of the earlier formed basin. Processes such as strain softening and mantle fertilization have a limited impact on the primary morphology and thermal state of rift systems before the crust is thinned to 10 km but may locally amplify relief and thermal heterogeneities.

SUMMARY Numerous processes such as metamorphic reactions, fluid and melt transfer and earthquakes occur at a subducting zone, but are still incompletely understood. These processes are affected, or even controlled, by the magnitude and distribution of stress and deformation mechanism. To eventually understand subduction zone processes, we quantify here stresses and deformation mechanisms in and around a subducting lithosphere, surrounded by asthenosphere and overlain by an overriding plate. We use 2-D thermomechanical numerical simulations based on the finite difference and marker-in-cell method and consider a 3200 km wide and 660 km deep numerical domain with a resolution of 1 km by 1 km. We apply a combined visco-elasto-plastic deformation behaviour using a linear combination of diffusion creep, dislocation creep and Peierls creep for the viscous deformation. We consider two end-member subduction scenarios: forced and free subduction. In the forced scenario, horizontal velocities are applied to the lateral boundaries of the plates during the entire simulation. In the free scenario, we set the horizontal boundary velocities to zero once the subducted slab is long enough to generate a slab pull force large enough to maintain subduction without horizontal boundary velocities. A slab pull of at least 1.8 TN m–1 is required to continue subduction in the free scenario. We also quantify along-profile variations of gravitational potential energy (GPE). We evaluate the contributions of topography and density variations to GPE variations across a subduction system. The GPE variations indicate large-scale horizontal compressive forces around the trench region and extension forces on both sides of the trench region. Corresponding vertically averaged differential stresses are between 120 and 170 MPa. Furthermore, we calculate the distribution of the dominant deformation mechanisms. Elastoplastic deformation is the dominant mechanism in the upper region of the lithosphere and subducting slab (from ca. 5 to 60 km depth from the top of the slab). Viscous deformation dominates in the lower region of the lithosphere and in the asthenosphere. Considering elasticity in the calculations has an important impact on the magnitude and distribution of deviatoric stress; hence, simulations with increased shear modulus, in order to reduce elasticity, exhibit considerably different stress fields. Limiting absolute stress magnitudes by decreasing the internal friction angle causes slab detachment so that slab pull cannot be transmitted anymore to the horizontal lithosphere. Applying different boundary conditions shows that forced subduction simulations are stronger affected by the applied boundary conditions than free subduction simulations. We also compare our modelled topography and gravity anomaly with natural data of seafloor bathymetry and free-air gravity anomalies across the Mariana trench. Elasticity and deviatoric stress magnitudes of several hundreds of MPa are required to best fit the natural data. This agreement suggests that the modelled flexural behaviour and density field are compatible with natural data. Moreover, we discuss potential applications of our results to the depth of faulting in a subducting plate and to the generation of petit-spot volcanoes.

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.

International audience; Constitutive models for rocks and soils typically incorporate some form of strain softening. Moreover, many plasticity models for frictional materials use a non-associated flow rule. Strain softening and non-associated flow rules can cause loss of well-posedness of the initial-value problem, which can lead to a severe mesh dependence in simulations and poor convergence of the iterative solution procedure. The inclusion of viscosity, which is a common property of materials, seems a natural way to restore well-posedness, but the mathematical properties of a rate-dependent model, and therefore the effectiveness with respect to the removal of mesh dependence , can depend strongly on how the viscous element is incorporated. Herein, we show that rate-dependent models which are commonly applied to problems in the Earth's lithosphere, such as plate tectonics, is very different from the approach typically adopted for more shallow geotechnical engineering problems. We analyse the properties of these models under dynamic loadings, using dispersion analyses and one-dimensional finite difference analyses, and complement them with two-dimensional simulations of a typical strain localisation problem under quasi-static loading conditions. Finally, we point out that a combined model, which features two viscous elements, may be the best way forward for modelling time-dependent failure processes in the deeper layers of the Earth, since it enables modelling of the creep characteristics typical of long-term behaviour, but also regularises the initial/boundary-value problem.

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.

International audience; Strain localization is a fundamental characteristic of plate tectonics. The resulting deformation structures shape the margins of continents and the internal structure of tectonic plates. To model the occurrence of faulting, geodynamic models generally rely on frictional plasticity. Frictional plasticity is normally embedded in visco‐plastic (V‐P) or visco‐elasto‐plastic (V‐E‐P) rheologies. This poses some fundamental issues, such as the difficulty, or often inability, to obtain a converged equilibrium state, and a severe grid sensitivity. Here, we study shear banding at crustal‐scale using a visco‐elasto‐viscoplastic (V‐E‐VP) model. We show that this rheology allows to accurately satisfy equilibrium, leads to shear band patterns that converge upon mesh refinement and preserves characteristic shear band angles. Moreover a comparison with analytic models and laboratory data reveals that V‐E‐VP rheology captures first order characteristics of frictional plasticity. V‐E‐VP models thus overcomes limitations of V‐P and V‐E‐P models and appears as an attractive alternative for geodynamic modeling.

International audience; The thermal state of the lithosphere is often considered as the main factor controlling the distribution of the structural and metamorphic features in tectonic systems (ancient versus modern tectonics). Deformation in ancient (weak and hot) lithospheres is distributed whereas deformation in modern (strong and cold) lithospheres is localised into prominent shear zones. The distributed deformation during the Precambrian suggests lower compressive strain rates are required, implying a long compression duration. However, the effects of variable strain rate on the compression of ancient lithospheres remain under-explored. A series of numerical models mimicking Precambrian lithospheric conditions is proposed and a broad range of thermal profiles and shortening rates is tested to investigate their influence on the resultant deformation modes. Two broad deformation styles were found to be significantly influenced by the magnitude of the shortening rates and stand out from the parametric study: (i) pop-downs and upper crustal thrusting at high strain rates that favour the formation of shear zones, and (ii) sedimentary cusping deformation at slower strain rates that generates well-defined vertical finite strain patterns. These two deformation modes were compared with natural examples, providing further insight into the evolution of their respective deformations. The vertical structures at the Hearne Craton (western Canadian Shield) could be explained via the sedimentary cusping mechanism by constraining our models with available geological data. Meanwhile, the vertical finite strain patterns from the proposed reference model now provide a tectonic template for future fluid-thermal convection simulations, to predict fluid flow and mineralisation potential.

Reliable numerical models of lithospheric deformation require robust solution methods. The latter should account for a complex and realistic rheological model and should also provide convergent and reproducible results.Here we present models of crustal-scale deformation that accurately capture the phenomenon of strain localisation in two-dimensions. The use of viscous regularisation yields convergent numerical results. We will compare linearisation methods (consistent tangent, effective viscosity) and discuss the implementation of rheological models (power-law viscous, hardening/softening laws). We will also present three-dimensional models of crustal-scale strain localisation that benefit from both the above-described methods and the computing power of graphical processing units (GPUs).

SUMMARY We present 2-D numerical simulations of convergence at a hyperextended passive margin with exhumed subcontinental mantle. We consider viscoelasto-plastic deformation, heat transfer and thermomechanical coupling by shear heating and associated thermal softening due to temperature dependent viscosity. The simulations show subduction initiation for convergence velocities of 2 cm yr−1, initial Moho temperatures of 525 °C and maximal deviatoric stresses of ca. 800 MPa, around the Moho, prior to localization. Subduction initiates in the region with thinned continental crust and is controlled by a thermally activated ductile shear zone in the mantle lithosphere. The shear zone temperature can be predicted with a recently published analytical expression. The criterion for subduction initiation is a temperature difference of at least 225 °C between predicted temperature and initial Moho temperature. The modelled forced subduction broadly agrees with geological data and reconstructions of subduction during closure of the Piemont-Liguria basin, caused by convergence of the European and Adriatic plates during the Alpine orogeny.

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.

International audience; Uranium mineralization is common in geological environments where oxidized uranium bearing fluids interact with reduced lithologies. This was the case for Variscan detachments in the South Armorican Massif, France. Numerous syn-tectonic peraluminous granites were emplaced in the footwall of detachments during the Carboniferous. Surface-derived fluids infiltrated the detachment planes and leached out the uranium from the magmatic uraninite that crystallized in these granites at the end of magmatic differentiation. On their way back to the surface, these fluids interacted with reducing black schists, which triggered precipitation of hydrothermal uraninite. This conceptual model serves as a basis for thermal hydraulic-mechanical modeling.

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; We propose an explanation to the enigmatic synrift erosional unconformities reported along the distal domain of several magma-poor rifted margins. Using thermomechanical numerical modeling, we show that transient emersion of (future) distal domains following a phase of subsidence can be explained by asynchronous necking of first the upper mantle and subsequently the crust, without the need of prominent normal faulting caused by strain softening, mantle phase transitions, or magmatic processes. When the upper crust and upper mantle are mechanically decoupled by a weak lower crust and, in the absence of any prominent rheological heterogeneity, upper mantle, necking starts first because of the higher deviatoric stresses associated with its larger effective viscosity. Consequently, the ductile lower crustal material flows toward the necked mantle domain, delaying thinning of the overlying crust. Once the necked lithospheric mantle has locally lost most of its strength, the overdeepened Moho moves upward toward an isostatically compensated depth. This flexural rebound causes uplift and emersion of distal parts of the rift system that are composed of still relatively thick crust and triggers the necking of the overlying crust. Early necking of the upper mantle causes a transient heating event with temperatures up to 750 • C at the base of the crust in the (future) distal domain. The onset of this thermal event slightly predates emersion of the (future) distal domain. These results are consistent with field observations and thermochronological data from the fossil Alpine Tethys margins, as well as with seismic observations from several present-day rifted margins.

International audience; The mechanisms responsible for the formation of extremely thinned continental crust (