Thibault Duretz

Discretization errors and free surface stabilization in the finite difference and marker‐in‐cell method for applied geodynamics: A numerical study

The finite difference–marker-in-cell (FD-MIC) method is a popular method in thermomechanical modeling in geodynamics. Although no systematic study has investigated the numerical properties of the method, numerous applications have shown its robustness and flexibility for the study of large viscous deformations. The model setups used in geodynamics often involve large smooth variations of viscosity (e.g., temperature-dependent viscosity) as well large discontinuous variations in material properties (e.g., material interfaces). Establishing the numerical properties of the FD-MIC and showing that the scheme is convergent adds relevance to the applications studies that employ this method. In this study, we numerically investigate the discretization errors and order of accuracy of the velocity and pressure solution obtained from the FD-MIC scheme using two-dimensional analytic solutions. We show that, depending on which type of boundary condition is used, the FD-MIC scheme is a second-order accurate in space as long as the viscosity field is constant or smooth (i.e., continuous). With the introduction of a discontinuous viscosity field characterized by a viscosity jump (η*) within the control volume, the scheme becomes first-order accurate. We observed that the transition from second-order to first-order accuracy will occur with only a small increase in the viscosity contrast (η* ≈ 5). We have employed two methods for projecting the material properties from the Lagrangian markers onto the Eulerian nodes. The methods are based on the size of the interpolation volume (4-cell, 1-cell). The use of a more local interpolation scheme (1-cell) decreases the absolute velocity and pressure discretization errors. We also introduce a stabilization algorithm that damps the potential oscillations that may arise from quasi free surface calculations in numerical codes that employ the strong form of the Stokes equations. This correction term is of particular interest for topographic modeling, since the surface of the Earth is generally represented by a free surface. Including the stabilization enables physically meaningful solutions to be obtained from our simulations, even in cases where the time step value exceeds the isostatic relaxation time. We show that including the stabilization algorithm in our FD stencil does not affect the convergence properties of our scheme. In order to verify our approach, we performed time-dependent simulations of free surface Rayleigh-Taylor instability.

Slab detachment in laterally varying subduction zones: 3-D numerical modeling

Understanding the three-dimensional (3-D) dynamics of subduction-collision systems is a longstanding challenge in geodynamics. We investigate the impact of slab detachment in collision systems that are subjected to along-trench variations. High-resolution thermomechanical numerical models, encompassing experimentally derived flow laws and a pseudo free surface, are employed to unravel lithospheric and topographic evolutions. First, we consider coeval subduction of adjacent continental and oceanic lithospheres (SCO). This configuration yields to two-stage slab detachment during collision, topographic buildup and extrusion, variable along-trench convergence rates, and associated trench deformation. The second setting considers a convergent margin, which is laterally limited by a transform boundary (STB). Such collisional system is affected by a single slab detachment, little trench deformation, and moderately confined upper plate topography. The effect of initial thermal slab age on SCO and STB models are explored. Similarities with natural analogs along the Arabia-Eurasia collision are discussed.

Shear heating-induced strain localization across the scales

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 level in the material outside the shear zone (stress drop of ). Conversely, Newtonian materials produce broad shear zones ( larger than for ) and exhibit a stronger thermal weakening (stress drops of ). Finally, we evaluate the thickness of shear zones for a wide range of strain rates (from to ) using both numerical models and an analytical scaling law. Our results suggest that thermally activated shear localization may occur from the scale of kilometre down to the nanometre.

From symmetric necking to localized asymmetric shearing: The role of mechanical layering

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.

The importance of structural softening for the evolution and architecture of passive margins

Lithospheric extension can generate passive margins that bound oceans worldwide. Detailed geological and geophysical studies in present and fossil passive margins have highlighted the complexity of their architecture and their multi-stage deformation history. Previous modeling studies have shown the significant impact of coarse mechanical layering of the lithosphere (2 to 4 layer crust and mantle) on passive margin formation. We built upon these studies and design high-resolution (~100–300 m) thermo-mechanical numerical models that incorporate finer mechanical layering (kilometer scale) mimicking tectonically inherited heterogeneities. During lithospheric extension a variety of extensional structures arises naturally due to (1) structural softening caused by necking of mechanically strong layers and (2) the establishment of a network of weak layers across the deforming multi-layered lithosphere. We argue that structural softening in a multi-layered lithosphere is the main cause for the observed multi-stage evolution and architecture of magma-poor passive margins.

M2Di: Concise and efficient MATLAB 2-D Stokes solvers using the Finite Difference Method

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.

Garnet xenocryst from petit-spot lavas as an indicator for off-axis mantle refertilization at intermediate spreading ridges

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