Dynamic modeling of sequences of earthquakes and aseismic slip (SEAS) provides a self-consistent, physics-based framework to connect, interpret, and predict diverse geophysical observations across spatial and temporal scales. Amid growing applications of SEAS models, numerical code verification is essential to ensure reliable simulation results but is often infeasible due to the lack of analytical solutions. Here, we develop two benchmarks for three-dimensional (3D) SEAS problems to compare and verify numerical codes based on boundary-element, finite-element, and finite-difference methods, in a community initiative. Our benchmarks consider a planar vertical strike-slip fault obeying a rate- and state-dependent friction law, in a 3D homogeneous, linear elastic whole-space or half-space, where spontaneous earthquakes and slow slip arise due to tectonic-like loading. We use a suite of quasi-dynamic simulations from 10 modeling groups to assess the agreement during all phases of multiple seismic cycles. We find excellent quantitative agreement among simulated outputs for sufficiently large model domains and small grid spacings. However, discrepancies in rupture fronts of the initial event are influenced by the free surface and various computational factors. The recurrence intervals and nucleation phase of later earthquakes are particularly sensitive to numerical resolution and domain-size-dependent loading. Despite such variability, key properties of individual earthquakes, including rupture style, duration, total slip, peak slip rate, and stress drop, are comparable among even marginally resolved simulations. Our benchmark efforts offer a community-based example to improve numerical simulations and reveal sensitivities of model observables, which are important for advancing SEAS models to better understand earthquake system dynamics.

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.