Numerical simulations of seismic cycles with rate-, state-, and temperature-dependent friction explain the various modes of seismic and aseismic ruptures in the brittle section of the lithosphere. However, the effects of viscoelastic flow in the ductile layers remain challenging to incorporate due to the wide range of length scales involved, from extremely localized within fault zones to widely distributed in the lower crust and asthenosphere. Here, we describe simulations of seismic cycles in a viscoelastic half-space using the integral method that combines discrete surface and volume elements to capture the coupling between brittle and ductile deformation. Viscoelastic flow is captured by cuboidal and tetrahedral volume elements within rectilinear and curvilinear meshes, respectively. The model resolves all phases of the seismic cycle under the radiation-damping approximation, including the nucleation and propagation of earthquake ruptures, but also the viscoelastic relaxation that follows in the ductile layers. We illustrate the approach in three dimensions with numerical simulations of seismic cycles on finite strike-slip and thrust faults overlying a viscoelastic lower crust with linear and nonlinear rheology. In two-dimensional models of subduction zones with the in-plane strain approximation, the ductile regions are meshed with triangle volume elements. The use of Green's functions only requires the discretization of the actively deforming region, resulting in a relatively small mesh. We provide open-source software implementing the method with parallel computing in a distributed architecture. The approach allows increasingly realistic representations of the lithosphere-asthenosphere system with nonlinear constitutive laws in structurally complex tectonic settings.

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.

The continental collision at the western boundary of the Indian continent formed the tectonically complex transpressional zones of the Sulaiman Fold Thrust (SFT) and Kirthar Fold Thrust (KFT) belts. Seismic hazard around the SFT is considered elevated, but shortening across its eastern boundary is poorly understood because of the scarcity of moderate-sized earthquakes in the last few decades. Here, we use Sentinel-1A ascending and descending interferometry to analyze the coseismic crustal deformation associated with the 2015 moment magnitude (Mw) 5.7 Dajal earthquake that occurred on the boundary thrust in the SFT belt. The surface displacement was caused by slip on a blind thrust and coseismic folding in the hanging wall. We use kinematic inversions to determine the distribution of slip on the frontal ramp and flexural slip along active axial surfaces for two end-member models of fault geometry. We first consider a double fault-bend fold system involving two sub-horizontal décollements separated by a ramp. Second, we consider a fault-propagation fold system where the frontal ramp terminates below a thick sediment layer. The fault-bend fold model includes slip on the décollement-ramp-décollement and flexural slip on two active axial surfaces initiated at the fault bends. For the fault-propagation fold, the model includes slip on the décollement-ramp system and flexural slip on the lower axial surface. In a preliminary step, the geometry of the ramp is optimized using a Monte Carlo method using a single asperity model. The geometry of the décollement and active axial surfaces is inferred based on balanced cross-sections, whereby the fold axes bisect the sediment layers across a fault bend. In both end-member models, a shallow décollement branches into a shallow ramp at approximately 7 km depth. However, the undeformed sediment of the overlying floodplain indicates that a ramp is the recent geomorphic feature. We conclude that the Dajal earthquake propagated along the base of the ramp, representing coseismic ramp failure over fault-propagation folding.