Finite thickness of shear bands in frictional viscoplasticity and implications for lithosphere dynamics

Permanent deformations in the lithosphere can occur in the brittle as well as in the ductile domain. For this reason, the inclusion of viscous creep and frictional plastic deformation is essential for geodynamic models. However, most currently available models of frictional plasticity are rate independent and therefore do not incorporate an internal length scale, which is an indispensible element for imposing a finite width of localized shear zones. Therefore, in computations of localization, either analytical or numerical, resulting shear zone widths tend to zero. In numerical computations, this manifests itself in a severe mesh sensitivity. Moreover, convergence of the global iterative procedure to solve the nonlinear processes is adversely affected, which negatively affects the reliability and the quality of predictions. The viscosity that is inherent in deformation processes in the lithosphere can, in principle, remedy this mesh sensitivity. However, elasto‐viscoplastic models that are commonly used in geodynamics assume a series arrangement of rheological elements (Maxwell‐type approach), which does not introduce an internal length scale. Here, we confirm that a different rheological arrangement that puts a damper in parallel to the plastic slider (Kelvin‐type approach) introduces an internal length scale. As a result, pressure and strain and strain rate profiles across the shear bands converge to finite values upon decreasing the grid spacing. We demonstrate that this holds for nonassociated plasticity with constant frictional properties and with material softening with respect to cohesion. Finally, the introduction of Kelvin‐type viscoplasticity also significantly improves the global convergence of nonlinear solvers

Toward robust and predictive geodynamic modeling : the way forward in frictional plasticity

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

Anisotropy and XKS-splitting from geodynamic models of double subduction: Testing the limits of interpretation

In this study, we develop three-dimensional geodynamic models to predict XKS-splitting for double subduction scenarios characterized by two outward dipping slabs. These models are highly relevant in various realistic settings, such as the central Mediterranean. We focus on the analysis of XKS-splitting, a key geophysical observable used to infer seismic anisotropy and mantle flow patterns predicted from these geodynamic models. Our geodynamic models simulate the concurrent subduction of two identical oceanic plates which are separated by a continental plate. The variation of the separating plate strength, cause a transition from a retreating to a stationary trench. The models provide detailed insights into the temporal evolution of mantle flow patterns, especially the amount of trench parallel flow, induced by these double subduction scenarios. In a second step, we use the well-known D-Rex model (Kaminski et al., 2004) to efficiently estimate the CPO development in response to plastic deformation produced by mantle flow. Based on the results of the D-Rex model, which includes the full elastic tensor of a deformed multiphase polycrystalline mantle aggregate within the three-dimensional model, we obtain synthetic apparent splitting parameters at receivers placed at the surface by applying multiple-layer anisotropic waveform modeling. Employing analytical techniques, we show the ambiguous nature of apparent splitting parameters, as already suggested by previous studies based on numerical modeling. In the light of the results, we postulate that a meaningful inversion, based on the commonly applied 2-layer anisotropic model, requires additional constraints on fast-axis orientation or strength of anisotropy (delay time). Finally, we show that constraints from our texture simulations (i.e., the integrated delay time) can be used to achieve unique 2-layer models that perfectly fit the synthetic observables. Such models could serve as reference for the interpretation of the observations. Our study highlights the necessity of combining geodynamic modeling and XKS-splitting analysis to shed light on complex upper mantle flow patterns such as those that might occur around subduction zones

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

A comparison of plasticity regularization approaches for geodynamic modeling

The emergence, geometry and activation of faults are intrinsically linked to frictional rheology. The latter is thus a central element in geodynamic simulations which aim at modeling the generation and evolution of fault zones and plate boundaries. However, resolving frictional strain localization in geodynamic models is problematic. In simulations, equilibrium cannot always be attained and results can depend on mesh resolution. Spatial and temporal regularization techniques have been developed to alleviate these issues. Herein, we investigate three popular regularization techniques, namely viscoplasticity, gradient plasticity and the use of a Cosserat continuum. These techniques have been implemented in a single framework based on an accelerated pseudo-transient solution strategy. The latter allows to explore the effects of regularization on shear banding using the same code and model configuration. We have used model configurations that involve three levels of complexity: from the emergence of a single isolated shear band to the visco-elasto-plastic stress buildup of a crust. All considered approaches allow to resolve shear banding, provide convergence upon mesh refinement and satisfaction of equilibrium. Viscoplastic regularization is straightforward to implement in geodynamic codes. Nevertheless, more stable shear banding patterns and strength estimates are achieved with computationally more expensive gradient and Cosserat-type regularizations. We discuss the relative benefits of these techniques and their combinations for geodynamic modeling. Emphasis is put on the potential of Cosserat-type media for geodynamic applications

Thinning mechanisms of heterogeneous continental lithosphere

The mechanisms responsible for the formation of extremely thinned continental crust (<10 km thick) and lithosphere during rifting remains debated. Observations from present-day and fossil passive margins highlight the role of deep-seated deformation, likely controlled by heterogeneities within the continental lithosphere, such as changing lithologies, mechanical anisotropies and inherited structures. We investigate the mechanisms of lithospheric thinning by exploring the role of pre-existing heterogeneities on the architecture and evolution of rifted margins. We estimate pre-rift pressure conditions (P0) vs. depth diagrams of crustal to lithospheric sections, to quantify rift-related modifications on inherited lithostatic pressure gradients. Two field examples from the Alpine Tethys margins in the Eastern and Southern Alps (SE Switzerland and N Italy) were selected to characterize: (1) the pre-rift architecture of the continental lithosphere; (2) the localization of rift-related deformation in distinct portions of the lithosphere; and (3) the interaction between pre-existing heterogeneities of the lithosphere and rift-related structures. These observations are compared with high-resolution, two-dimensional thermo-mechanical numerical models. The design of the models takes into account pre-existing mechanical heterogeneities representing the initial pre-rift architecture of the continental lithosphere. Extensional structures consist of high-angle and low-angle normal faults, anastomosing shear-zones and decoupling horizons. Such structures accommodate the lateral extraction of mechanically stronger levels derived from the middle to lower crust. As a result, the extremely thinned continental crust in Tethyan passive margins represents the juxtaposition and amalgamation of distinct strong levels of the crust separated by major extensional structures identified by sharp pressure gradients. Future work should determine the applicability of these results to other present-day and fossil rifted margins

A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air' method

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

Magma differentiation in dykes: from field evidence to numerical study

Mechanisms active during magma ascent in dykes are not still well quantified. From field observations, it is apparent that most dykes actually contain a crystalline load. The presence of a crystalline load modifies the effective rheology of such a system and thus the flow behaviour. For example, increasing both the density and viscosity of each crystal, compared to the melt, will cause a reduction of the ascent velocity and will modify the shape of the velocity profile from “Poiseuille like” to “Bingham-like”. One issue related to the introduction of a crystalline load concerns the possibility for crystals to be segregated from a viscous granitic melt phase during magma ascent. The implications of such a process on magmatic differentiation have not previously been considered, nor has such a process been previously investigated via numerical models. In this study, we examine the flow dynamics of a poly-crystal bearing granitic melt ascending in a dyke via numerical models. Results showed that the melt phase can be squeezed out from a crystal-rich magma when subjected to a given pressure gradient range and that crystals behaviour is strongly dependant on their size and density. This demonstrates that crystal-melt segregation in dykes during granitic magma ascent constitutes a viable mechanism for magmatic differentiation. In order to quantify such a mechanism, results have been compared to two well-characterized granites from the Armorican Massif (France). Geochemically, these two granites originate from the same magmatic source, but the shallowest one (Questembert granite) can be derived from the deepest (Lizio granite) by a larger amount of crystal segregation during magma ascent in dykes. In our numerical models, we therefore use as initial conditions, the volume fraction and geochemical properties of the “Lizio-type” magma and compute the vertical distance required to end up with the “Questembert-type” magma. The comparison of the results allowed us to estimates the parameters such as pressure gradient, distance, viscosity, etc. . . needed, to explain their different degrees of geochemical differentiation

Numerical modelling of magma transport in dykes

The rheology and dynamics of an ascending pure melt in a dyke have been extensively studied in the past. From field observations, it is apparent that most dykes actually contain a crystalline load. The presence of a crystalline load modifies the effective rheology of such a system and thus the flow behaviour. Indeed, the higher density and viscosity of each crystal, compared to the melt, cause a decrease of the ascent velocity and modify the shape of the velocity profile, from a typical Poiseuille flow, to a Bingham-type flow. A common feature observed in the field is the arrangement of crystals parallel or at a very low angle to the edge of the dyke. Such a structural arrangement is often interpreted as the result of magma flow, which caused the crystals to rotate and align within the flow direction, but this process remains unclear. Another issue related to the introduction of a crystalline load concerns the possibility for crystals to be segregated from a viscous granitic melt phase during magma ascent. The implications of such a process on magmatic differentiation have not previously been considered, nor has such a process been previously investigated via numerical models. In this study, we examine the flow dynamics of a crystal bearing granitic melt ascending in a dyke via numerical models. In our models, both the crystal and melt phases are represented as highly viscous fluids in a Stokes regime. Our results reveal that the presence of crystals in the melt modifies the magma velocity profile across the dyke. Furthermore, we observe that whilst crystals continually rotate in the shear flow, over one period of revolution, their major axis has a high probability to be aligned parallel to the flow direction. Moreover, some experiments showed that the melt phase can effectively be squeezed out from a crystal-rich magma when subjected to a given pressure gradient range. This demonstrates that crystal-melt segregation in dykes during granitic magma ascent constitutes a viable mechanism for magmatic differentiation