A rapid numerical routine for viscoelastic earthquake cycle simulations.

Earth’s largest quakes and trans-oceanic tsunamis emanate from subduction zones around the world. Following such large earthquakes, viscoelastic processes and on-fault aseismic fault slip play a crucial role in dissipating the stresses induced by the earthquake, facilitating the solid Earth's return to equilibrium.  The rheological properties of lithospheric rocks govern these postseismic processes and influence time-dependent deformation during the earthquake cycle. Geodetic observations offer an opportunity to constrain these rheological properties, providing valuable insights into the regional lithospheric structure, and potentially improving our understanding of earthquake-related hazards.   

To build intuition for geodetically recorded postseismic deformation, we develop a robust and efficient two-dimensional quasi-static periodic earthquake cycle simulator exploiting the boundary element method and semi-analytical solutions to systems of coupled ordinary differential equations. We investigate the impact of lateral and depth-dependent variations in the viscosity structure of the mantle wedge and the oceanic mantle, to discern their respective contributions and roles in surface deformation observations. We account for the long-term viscous flow rate in the mantle and show that neglecting this term in the earthquake cycle introduces biases in the effective viscosity structure of the lithosphere-asthenosphere system, particularly in the context of power-law rheologies. The low computational cost of our numerical routine makes it ideal for incorporating into future inverse modelling frameworks to estimate regional rheological structure from geodetic observations of subduction zone earthquake cycles.