Stress and deformation mechanisms at a subduction zone: insights from 2-D thermomechanical numerical modelling

SUMMARY Numerous processes such as metamorphic reactions, fluid and melt transfer and earthquakes occur at a subducting zone, but are still incompletely understood. These processes are affected, or even controlled, by the magnitude and distribution of stress and deformation mechanism. To eventually understand subduction zone processes, we quantify here stresses and deformation mechanisms in and around a subducting lithosphere, surrounded by asthenosphere and overlain by an overriding plate. We use 2-D thermomechanical numerical simulations based on the finite difference and marker-in-cell method and consider a 3200 km wide and 660 km deep numerical domain with a resolution of 1 km by 1 km. We apply a combined visco-elasto-plastic deformation behaviour using a linear combination of diffusion creep, dislocation creep and Peierls creep for the viscous deformation. We consider two end-member subduction scenarios: forced and free subduction. In the forced scenario, horizontal velocities are applied to the lateral boundaries of the plates during the entire simulation. In the free scenario, we set the horizontal boundary velocities to zero once the subducted slab is long enough to generate a slab pull force large enough to maintain subduction without horizontal boundary velocities. A slab pull of at least 1.8 TN m–1 is required to continue subduction in the free scenario. We also quantify along-profile variations of gravitational potential energy (GPE). We evaluate the contributions of topography and density variations to GPE variations across a subduction system. The GPE variations indicate large-scale horizontal compressive forces around the trench region and extension forces on both sides of the trench region. Corresponding vertically averaged differential stresses are between 120 and 170 MPa. Furthermore, we calculate the distribution of the dominant deformation mechanisms. Elastoplastic deformation is the dominant mechanism in the upper region of the lithosphere and subducting slab (from ca. 5 to 60 km depth from the top of the slab). Viscous deformation dominates in the lower region of the lithosphere and in the asthenosphere. Considering elasticity in the calculations has an important impact on the magnitude and distribution of deviatoric stress; hence, simulations with increased shear modulus, in order to reduce elasticity, exhibit considerably different stress fields. Limiting absolute stress magnitudes by decreasing the internal friction angle causes slab detachment so that slab pull cannot be transmitted anymore to the horizontal lithosphere. Applying different boundary conditions shows that forced subduction simulations are stronger affected by the applied boundary conditions than free subduction simulations. We also compare our modelled topography and gravity anomaly with natural data of seafloor bathymetry and free-air gravity anomalies across the Mariana trench. Elasticity and deviatoric stress magnitudes of several hundreds of MPa are required to best fit the natural data. This agreement suggests that the modelled flexural behaviour and density field are compatible with natural data. Moreover, we discuss potential applications of our results to the depth of faulting in a subducting plate and to the generation of petit-spot volcanoes.