A fast numerical routine to model viscoelastic earthquake cycles in subduction zones

Abstract Understanding the mechanics and rheological properties of the solid Earth requires integrating observations of earthquake cycle deformation with advanced numerical models of the underlying processes. In this article, we introduce a new numerical modeling tool designed for two-dimensional subduction zone geometries that solves the equations for mechanical equilibrium in the lithosphere–asthenosphere system to predict displacements, strains, and stresses throughout the medium. Using an extended Boundary Element Method (BEM) framework, we reduce the governing partial differential equations to a set of coupled ordinary differential equations, simulating periodic earthquake cycles in a viscoelastic medium with spatially variable power-law viscous rheologies and rate-dependent fault friction. Our approach overcomes key limitations of existing BEM models by ensuring: (1) mechanical consistency across timescales, from individual earthquake cycles to geological periods, and (2) precise stress transfer calculations in a viscoelastic medium. Crucially, we demonstrate that while consistent long-term strain rates are essential for predicting postseismic deformation, a Newtonian corner-flow approximation of the background flow field provides a highly accurate and computationally efficient alternative to more computationally demanding non-linear solutions. We also show that for spatially heterogeneous linear viscoelastic materials, the coupled ordinary differential equations can be simplified to an eigenvalue problem, significantly enhancing computational efficiency. These advancements offer a powerful tool for predicting spatio-temporal patterns of surface displacement given complex mantle structures and lay the foundation for high-dimensional inverse problems to infer constitutive properties of the lithosphere–asthenosphere system.