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 a 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. We also demonstrate 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 spatiotemporal 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.

Surface deformation associated with continental earthquake ruptures includes localized deformation on the faults, as well as deformation in the surrounding medium through distributed and/or diffuse processes. However, the role of the diffuse part of the surface deformation to the overall rupture process, as well as its underlying physical mechanisms are not yet well understood. In this study, we compute high-resolution near-fault displacement maps from optical image correlations for the 2021/05/21 Mw7.4 Maduo, Tibet, strike-slip earthquake, and measure the contributions of the different deformation components to the surface deformations for that event. Results show that surface slip along primary faults accommodates, on average, only ~25% of the total surface deformation. Majority of the surface coseismic deformation is in fact accommodated by diffuse deformation,especially in the epicentral area where no surface slip was observed. In fact, the contribution of the diffuse deformation increases as localized deformation on the fault decreases. Localized deformation also decreases with decreasing total surface displacement. These observations highlight a gradual localization of the surface coseismic deformation, from regions of diffuse low (0.1-0.3%) strain, to regions of highly localized (>1 %) strain, with increasing coseismic displacement. Using simple two-dimensional mechanical models we show that diffuse deformation may correspond to elastoplastic bulk yielding, accounting for the deficit in shallow fault slip in the regions of surface rupture gap.

Understanding the nature and behavior of the rocks in boundary zones between tectonic plates is important to improve our understanding of earthquake-associated hazard. Laboratory experiments can derive models that explain material behavior on small scales and under controlled conditions. These models can also be tested on observations of surface motion near plate boundaries: Fitting surface displacements from earthquakes (either shortterm offsets or longterm motion) yields estimates of rock properties for each model. However, using only observations from a single earthquake (from immediately after the quake and/or the subsequent years), may not allows us to confidently distinguish between models. In this study, we investigate the potential of using the displacement timeseries from multiple earthquakes, as well as the period between the quakes, to distinguish between proposed models. We use methods that enable comparison between models and parameters taking into account uncertainties, and perform our assessment on an artificial dataset.

Viscoelastic processes in the upper mantle redistribute seismically generated stresses and modulate crustal deformation throughout the earthquake cycle. Geodetic observations of these motions at the Earth’s surface offer the possibility of constraining the rheology of the upper mantle. Parsimonious representations of viscoelastically modulated deformation should simultaneously be able to explain geodetic observations of rapid postseismic deformation and near-fault strain localization late in the earthquake cycle. We compare predictions from time-dependent forward models of deformation over the entire earthquake cycle on and surrounding an idealized vertical strike-slip fault in a homogeneous elastic crust underlain by a homogeneous viscoelastic upper mantle. We explore three different rheologies as inferred from laboratory experiments: 1) linear-Maxwell, 2) linear-Burgers, 3) power-law. Both the linear Burgers and power-law rheological models can be made consistent with fast and slow deformation phenomenology from across the entire earthquake cycle, while the single-layer linear Maxwell model cannot. The kinematic similarity of linear Burgers and power-law models suggests that geodetic observations alone are insufficient to distinguish between them, but indicate that one may serve as a proxy for the other. However, the power-law rheology model displays a postseismic response that is strongly earthquake magnitude dependent, which may offer a partial explanation for observations of limited postseismic deformation near magnitude 6.5-7.0 earthquakes. We discuss the role of mechanical coupling between frictional slip and viscous creep in controlling the time-dependence of regional stress transfer following large earthquakes and how this may affect the seismic hazard and risk to communities living close to fault networks.