Benchmark of three-dimensional numerical models of subduction against a laboratory experiment

In this paper, we validate three-dimensional numerical models of subduction against a laboratory experiment. The design of our numerical models follows that of the laboratory, and employs a single subducting plate fixed at its trailing edge without an overriding plate. Rheology and buoyancy are not thermally dependent, but importantly the model has a free surface boundary. All the models include a sticky air layer, which is an approximation to a free surface boundary condition, as opposed to a free-slip condition that unrealistically binds the plate to the upper surface. We use a parallel, Python, particle-in-cell, finite-element scheme, and present 8 different models that differ in their viscosity at the interfaces of the subducting plate. We consider harmonic and arithmetic averaging of the viscosity together with a plate composite rheology associated with arithmetic averaging, which similarly to the harmonic averaging results in two lubricated sub-layers at the plate lower and upper interfaces. Despite showing a similar subduction-induced poloidal and toroidal flow in the mantle, the agreement of the models with the experiment is overall ordinary. The rates at which the slab tip sinks and the plate hinge retreats are approaching the laboratory counterparts to within 7-18%, at best, when using the composite plate models. Errors in the characterization of the laboratory parameters can account for 5% of such differences. As the composite plate sub-layers are made stiffer, the plate hinge also begins to depart from the steady retreat observed in the laboratory experiment. All numerical models exhibit a secondary Rayleigh-Taylor instability along the plate edge that is not observed in the laboratory experiment. The development of this instability is further delayed when the numerical spatial resolution is increased, or when the composite plate sub-layers are made stiffer. This benchmark exercise shows that 1) the glucose surface tension in the laboratory, not taken into account in the numerical modelling, is sufficient to prevent the flank instability during the retreat, and 2) numerical models are usually free from it due to their use of stiffer plate sub-layers and/or their display of a limited amount of retreat compared to this study.