Evaluation of the presence of a weak layer in the numerical simulation of lithospheric subduction

One challenge to numerically simulate the subduction of cold oceanic lithosphere under continental lithosphere is the preservation of the decoupling between the subducting and upper plates for tens of millions of years. One strategy to simulate the persistence of the decoupling is the continuous entrainment of a weak layer (i.e. with low effective viscosity) at the top of the oceanic plate, representing a lubrication between both plates. However, variations on the thickness and rheological structure of this weak layer affect the geodynamic evolution of the subducting plate, modifying the geometry and degree of interactions between the lithospheric plates. 

In the present work we evaluated how the variation of the geometry, viscosity and density of the weak layer, relative to the surrounding lithosphere, can affect the lubrication between the two lithospheric plates. We performed a series of 2D numerical simulations using a finite element code for thermochemical convection. The code solves the Stokes flow for a fluid using the Boussinesq approximation in a Cartesian coordinate system, considering that the viscosity varies exponentially as a function of the temperature. In the present visco-plastic approach, the effective viscosity is determined by the combined effect of a viscous component, assuming the Frank-Kamenetskii rheology, and plastic deformation, following the Byerlee's friction law. 

In our numerical scenarios, the subduction is produced by the negative buoyancy of the cold oceanic lithosphere, without the imposition of an external velocity as boundary conditions. The time range of the simulation is of the order of 50 million years. In the initial simulation, a weak zone is imposed in the region between the two plates. This zone presents low viscosity and density relative to the surrounding lithosphere. As the oceanic slab is subducted, the weak zone is deformed and dragged. This removes the lubrication until utterly coupling the lithospheric plates, generating the thickening of the continental lithosphere below the trench region. To preserve the decoupling along all the simulation time, an extra continuous weak layer on top of the oceanic plate is added with low density and viscosity. In this scenario, the first weak zone is still dragged by the subducting plate, but the additional weak layer keeps a lubrication zone between the plates, preventing the coupling of the two lithospheric plates. Therefore, adding a continuous weak layer on top of the oceanic crust together with a weak zone prevents the coupling of the subducting and overriding plates when the effective viscosity of the weak layer is smaller than ~10¹⁹ Pa s. These numerical scenarios are used to analyse the subduction pattern of the Nazca plate observed in the southeastern portion of South America, using as constraints the slab geometry of the subducting oceanic plate derived from the Slab2 model.