A benchmark comparison of spontaneous subduction models – towards a free surface

Numerically modelling the dynamics of a self-consistently subducting lithosphere is a challenging task because of the decoupling problems of the slab from the free surface. We address this problem with a benchmark comparison between various numerical codes (Eulerian and Lagrangian, Finite Element and Finite Difference, with and without markers) as well as a laboratory experiment. The benchmark test consists of three prescribed setups of viscous flow, driven by compositional buoyancy, and with a low viscosity, zero-density top layer to approximate a free surface. Alternatively, a fully free surface is assumed. Our results with a weak top layer indicate that the convergence of the subduction behaviour with increasing resolution strongly depends on the averaging scheme for viscosity near moving rheological boundaries. Harmonic means result in fastest subduction, arithmetic means produces slow subduction and geometric mean results in intermediate behaviour. A few cases with the infinite norm scheme have been tested and result in convergence behaviour between that of arithmetic and geometric averaging. Satisfactory convergence of results is only reached in one case with a very strong slab, while for the other cases complete convergence appears mostly beyond presently feasible grid resolution. Analysing the behaviour of the weak zero-density top layer reveals that this problem is caused by the entrainment of the weak material into a lubrication layer on top of the subducting slab whose thickness turns out to be smaller than even the finest grid resolution. Agreement between the free surface runs and the weak top layer models is satisfactory only if both approaches use high resolution. Comparison of numerical models with a free surface laboratory experiment shows that (1) Lagrangian-based free surface numerical models can closely reproduce the laboratory experiments provided that sufficient numerical resolution is employed and (2) Eulerian-based codes with a weak surface layer reproduce the experiment if harmonic or geometric averaging of viscosity is used. The harmonic mean is also preferred if circular high viscosity bodies with or without a lubrication layer are considered. We conclude that modelling the free surface of subduction by a weak zero-density layer gives good results for highest resolutions, but otherwise care has to be taken in (1) handling the associated entrainment and formation of a lubrication layer and (2) choosing the appropriate averaging scheme for viscosity at rheological boundaries.

[1]  G. Marquart,et al.  The influence of second-scale convection on the thickness of continental lithosphere and crust , 1991 .

[2]  Louis Moresi,et al.  Numerical investigation of 2D convection with extremely large viscosity variations , 1995 .

[3]  R. W. Griffiths,et al.  Laboratory models of the thermal evolution of the mantle during rollback subduction , 2003, Nature.

[4]  W. Schellart.  Quantifying the net slab pull force as a driving mechanism for plate tectonics , 2004 .

[5]  H. Schmeling,et al.  The influence of olivine metastability on deep subduction of oceanic lithosphere , 2000 .

[6]  W. Schellart.  Influence of the subducting plate velocity on the geometry of the slab and migration of the subduction hinge , 2005 .

[7]  Roberto F. Weinberg,et al.  Polydiapirs: multiwavelength gravity structures , 1992 .

[8]  P. Raviart.  Finite element methods and Navier-Stokes equations , 1979 .

[9]  David A. Yuen,et al.  Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties , 2003 .

[10]  M. Moroni,et al.  Mapping mantle flow during retreating subduction: Laboratory models analyzed by feature tracking , 2006 .

[11]  Francesca Funiciello,et al.  Role of lateral mantle flow in the evolution of subduction systems: insights from laboratory experiments , 2004 .

[12]  S. Schmalholz,et al.  Numerical simulations of parasitic folding in multilayers , 2006 .

[13]  G. Morra,et al.  Dynamic models of downgoing plate-buoyancy driven subduction: subduction motions and energy dissipation , 2007 .

[14]  Boris J. P. Kaus,et al.  Comparison of Eulerian and Lagrangian numerical techniques for the Stokes equations in the presence of strongly varying viscosity , 2008 .

[15]  D. May,et al.  Evolution and diversity of subduction zones controlled by slab width , 2007, Nature.

[16]  Marcin Dabrowski,et al.  MILAMIN: MATLAB‐based finite element method solver for large problems , 2008 .

[17]  S. Labanieh,et al.  Dynamical effects of subducting ridges: insights from 3-D laboratory models , 2005, 1011.4041.

[18]  Louis-Noel Moresi,et al.  Influence of trench width on subduction hinge retreat rates in 3‐D models of slab rollback , 2006 .

[19]  Marie-Pierre Doin,et al.  A comparison of methods for the modeling of thermochemical convection , 1997 .

[20]  J. Z. Zhu,et al.  The finite element method , 1977 .

[21]  D. Giardini,et al.  The development of slabs in the upper mantle: Insights from numerical and laboratory experiments , 1999 .

[22]  P. Davy,et al.  Experiments on shortening of a 4-layer model of the continental lithosphere , 1991 .

[23]  B. Stöckhert,et al.  Exhumation of high‐pressure metamorphic rocks in a subduction channel: A numerical simulation , 2002 .

[24]  W. Schellart.  Kinematics of subduction and subduction-induced flow in the upper mantle , 2004 .

[25]  Boris J. P. Kaus,et al.  Effects of elasticity on the Rayleigh-Taylor instability: implications for large-scale geodynamics , 2005 .

[26]  D. Giardini,et al.  Dynamics of retreating slabs : 2. Insights from three-dimensional laboratory experiments , 2003 .

[27]  P. Olson,et al.  A laboratory model of subduction zone anisotropy , 1998 .

[28]  Gabriele Morra,et al.  Dynamics of retreating slabs: 1. Insights from two-dimensional numerical experiments , 2003 .

[29]  D. Rubie,et al.  The influence of olivine metastability on the dynamics of subduction , 1999 .

[30]  D. Giardini,et al.  Episodic Back-arc extension during restricted mantle convection in the Central Mediterranean , 2001 .

[31]  Dinesh Manocha,et al.  Applied Computational Geometry Towards Geometric Engineering , 1996, Lecture Notes in Computer Science.

[32]  H. Schmeling,et al.  On the effects of the lithosphere on mantle convection and evolution , 1982 .

[33]  David A. Yuen,et al.  Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems , 2007 .

[34]  R. Hilst,et al.  A laboratory investigation of effects of trench migration on the descent of subducted slabs , 1995 .

[35]  W. Jacoby.  Paraffin model experiment of plate tectonics , 1976 .

[36]  K. Bathe.  Finite Element Procedures , 1995 .

[37]  P Cundall,et al.  A MICROCOMPUTER PROGRAM FOR MODELLING LARGE-STRAIN PLASTICITY PROBLEMS. PROCEEDINGS OF THE SIXTH INTERNATIONAL CONFERENCE ON NUMERICAL METHODS IN GEOMECHANICS, 11-15 APRIL 1988, INNSBRUCK, AUSTRIA. VOLUMES 1 - 3 , 1988 .

[38]  Robert B. Trumbull,et al.  Numerical models of crustal scale convection and partial melting beneath the Altiplano–Puna plateau , 2002 .

[39]  P. Olson,et al.  Laboratory experiments on the structure of subducted lithosphere , 1995 .

[40]  Wolfgang R. Jacoby,et al.  On modelling the lithosphere in mantle convection with non-linear rheology , 1982 .

[41]  Ross W. Griffiths,et al.  Folding of Viscous Plumes Impinging On A Density Or Viscosity Interface , 1988 .

[42]  P. A. Cundall,et al.  An Explicit Inertial Method for the Simulation of Viscoelastic Flow: An Evaluation of Elastic Effects on Diapiric Flow in Two- and Three- Layers Models , 1993 .

[43]  S. Schmalholz.  Scaled amplification equation: A key to the folding history of buckled viscous single-layers , 2006 .

[44]  David A. Yuen,et al.  Rayleigh^Taylor instabilities from hydration and melting propel 'cold plumes' at subduction zones , 2003 .

[45]  R. Weijermars,et al.  Scaling of Newtonian and non-Newtonian fluid dynamics without inertia for quantitative modelling of rock flow due to gravity (including the concept of rheological similarity) , 1986 .

[46]  A. I. Shemenda,et al.  Physical modelling of lithosphere subduction in collision zones , 1992 .

[47]  M. Gurnis,et al.  Constraints on the lateral strength of slabs from three-dimensional dynamic flow models , 1996 .

[48]  D. Sulsky.  Erratum: Application of a particle-in-cell method to solid mechanics , 1995 .

[49]  Louis Moresi,et al.  Role of temperature‐dependent viscosity and surface plates in spherical shell models of mantle convection , 2000 .

[50]  Chris Kincaid,et al.  An experimental study of subduction and slab migration , 1987 .

[51]  N. Bellahsen,et al.  Dynamics of Subduction and Plate Motion in Laboratory Experiments. , 2004 .