Lagrangian Free Subduction Model Based on Viscoelastic Material using COMSOL Multiphysics®
Subduction, the process through which the cold and dense oceanic lithosphere sinks into the highly viscous mantle, is the main driving force of plate tectonics and is important for recycling the Earth's materials. Active subduction has not been observed in any other planet in our solar system; it is an indicator for classifying stagnant-lid and plate-like planets. For several decades, numerical modeling studies based on the finite element method have been performed to better understand the dynamics of subducting slabs. A 2D Lagrangian free subduction model using Maxwell viscoelastic materials was built primarily based on a commercial software called ABAQUS®. To improve the reliability of numerical models, benchmark studies should be performed with other software. For instance, the finite element package COMSOL Multiphysics®, which provides a user-friendly GUI and massive parallelization calculations, can be freely applied to various physical properties including geometry, element size, and shape function order. Studies looking at subduction models based on purely viscous fluids and governed by the Stokes equation in the Eulerian grid, have been conducted using COMSOL Multiphysics®; however, Lagrangian subduction modeling based on solid mechanics has never been quantitatively evaluated. Further, previous studies have assumed that the mechanical strengths of oceanic plates are simply iso- or layered viscosity. Contrary to this assumption, in order to model the realistic viscosity structure inside a slab, both nonlinear creep and brittle rheology must be applied, which are temperature and pressure dependent. In this study, we developed both 2D and 3D fully-coupled thermal (e.g. Heat Transfer interface), mechanical (e.g. Solid Mechanics interface), free subduction models that reflect realistic rheology and successfully simulate subduction self-consistently driven by gravity only. We also used the element type of free quad i.e. shape function of quadratic serendipity. Our models showed slab behavior and morphology evolving over time, with high stress concentrations in the upper and lower hinges of the slab. Further, we conducted numerical experiments on the temperature distribution in the slab during the subduction process using a Heat Transfer in Solid interface. We derived the cold slab core tongue structure by introducing the half space cooling model as the initial condition and the mantle geotherm as the Dirichlet boundary condition. Our model includes a realistic strength envelope curve for inside the slab obtained by applying creep rheology dependent on the temperature and stress, and derived from a fully coupled thermo-mechanical model. In subduction dynamics, the role of the Winkler foundation, used to mimic the isostatic restoring force caused by the low density of trench infill material (e.g., sediment and water), is still unclear. We quantified the effect of the Winker foundation due to its nature as an essential analytical element in the Lagrangian free subduction model. In other words, we suggest that if sedimentary accumulation rates vary gradually along trench axis, differential trench rollback rates could be obtained and oblique subduction could occur. Our results, including the geodynamic implication that the Earth’s surface process contributes to subduction physics, are expected to improve the understanding of the structure of deep subduction.