Improving subduction interface implementation in dynamic numerical models

Numerical subduction models often implement an entrained weak layer (WL) to facilitate decoupling of the slab and upper plate. This approach is attractive in its simplicity, and can provide stable, asymmetric subduction systems that persist for many tens of millions of years. In this study we undertake a methodological analysis of the WL approach, and use these insights to guide improvements to the implementation. The issue that primarily motivates the study is the emergence of significant spatial and temporal thickness variations within the WL. We show that these variations are mainly the response to volumetric flux gradients, caused by the change in boundary conditions as the WL material enters and exits the zone of decoupling. The time taken to reach a quasi-equilibrium thickness profile will depend on the total plate convergence, and is around 7 Myr for the models presented here. During the transient stage, width variations along the WL can exceed 4×, which may impact the effective strength of the interface, through physical effects if the rheology is linear, or simply if the interface becomes inadequately numerically resolved. The transient stage also induces strong sensitivity to model resolution. By prescribing a variable-thickness WL at the outset of the model, and by controlling the limits of the layer thickness during the model evolution, we find improved stability and resolution convergence of the models.