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 overwhelming observational difficulties and the complexity of earthquake physics have rendered seismic hazard assessment largely empirical. Despite increasingly high-quality geodetic, seismic, and field observations, data-driven earthquake imaging yields stark differences and physics-based models explaining all observed dynamic complexities are elusive. Here we present data-assimilated 3D dynamic rupture models which untwine California's biggest earthquakes in more than 20 years: the moment magnitude (M w ) 6.4 Searles Valley and M w 7.1 Ridgecrest, California, sequence breaking multiple segments of the same fault system. Our models use supercomputing to find the link between the two large earthquakes. We unify the uniquely high-quality strong-motion and teleseismic, field mapping, high-rate GNSS, and space geodetic foreshock and mainshock datasets with earthquake physics. We find that the regional structure, the ambient long-and short-term stress, as well as the dynamic and static fault system interactions, are conjointly crucial to understand the dynamics and delays of the sequence.Dynamic rupture of a statically strong yet dynamically weak fault system is driven by overpressurized fluids and low dynamic friction in our models. The observed earthquake complexity results from static and dynamic stress changes acting across a non-vertical quasi-orthogonal conjugate fault structure. We demonstrate that joint physics-based and data-driven illumination of the mechanics of complex fault systems and earthquake sequences is possible when reconciling dense earthquake recordings, 3D regional structure and stress models. We foresee that physics-based interpretation of big observational data-sets characterizing complex nonlinear systems will have a transformative impact on future geohazard mitigation.

Numerical modeling of earthquake dynamics and derived insight for seismic hazard relies on credible, reproducible model results. The sequences of earthquakes and aseismic slip (SEAS) initiative has set out to facilitate community code comparisons, and verify and advance the next generation of physics-based earthquake models that reproduce all phases of the seismic cycle. With the goal of advancing SEAS models to robustly incorporate physical and geometrical complexities, here we present code comparison results from two new benchmark problems: BP1-FD considers full elastodynamic effects, and BP3-QD considers dipping fault geometries. Seven and eight modeling groups participated in BP1-FD and BP3-QD, respectively, allowing us to explore these physical ingredients across multiple codes and better understand associated numerical considerations. With new comparison metrics, we find that numerical resolution and computational domain size are critical parameters to obtain matching results. Codes for BP1-FD implement different criteria for switching between quasi-static and dynamic solvers, which require tuning to obtain matching results. In BP3-QD, proper remote boundary conditions consistent with specified rigid body translation are required to obtain matching surface displacements. With these numerical and mathematical issues resolved, we obtain excellent quantitative agreements among codes in earthquake interevent times, event moments, and coseismic slip, with reasonable agreements made in peak slip rates and rupture arrival time. We find that including full inertial effects generates events with larger slip rates and rupture speeds compared to the quasi-dynamic counterpart. For BP3-QD, both dip angle and sense of motion (thrust versus normal faulting) alter ground motion on the hanging and foot walls, and influence event patterns, with some sequences exhibiting similar-size characteristic earthquakes, and others exhibiting different-size events. These findings underscore the importance of considering full elastodynamics and nonvertical dip angles in SEAS models, as both influence short- and long-term earthquake behavior and are relevant to seismic hazard.

The 2016 Central Italy earthquake sequence is characterized by remarkable rupture complexity, including highly heterogeneous slip across multiple faults in an extensional tectonic regime. The dense coverage and high quality of geodetic and seismic data allow to image intriguing details of the rupture kinematics of the largest earthquake of the sequence, the Mw 6.5 October 30th, 2016 Norcia earthquake, such as an energetically weak nucleation phase. Several kinematic models suggest multiple fault planes rupturing simultaneously, however, the mechanical viability of such models is not guaranteed.Using 3D dynamic rupture and seismic wave propagation simulations accounting for two fault planes, we constrain 'families' of spontaneous dynamic models informed by a high-resolution kinematic rupture model of the earthquake. These families differ in their parameterization of initial heterogeneous shear stress and strength in the framework of linear slip weakening friction.First, we dynamically validate the kinematically inferred two-fault geometry and rake inferences with models based on only depth-dependent stress and constant friction coefficients. Then, more complex models with spatially heterogeneous dynamic parameters allow us to retrieve slip distributions similar to the target kinematic model and yield good agreement with seismic and geodetic observations. We discuss the consistency of the assumed constant or heterogeneous static and dynamic friction coefficients with mechanical properties of rocks at 3-10 km depth characterizing the Italian Central Apennines and their local geological and lithological implications. We suggest that suites of well-fitting dynamic rupture models belonging to the same family generally exist and can be derived by exploiting the trade-offs between dynamic parameters.Our approach will be applicable to validate the viability of kinematic models and classify spontaneous dynamic rupture scenarios that match seismic and geodetic observations at the same time as geological constraints.

We design the first fully dynamic 3-D earthquake simulations based on geodetic coupling models for the Cascadia megathrust. • Segmentation in the stress drop is needed to produce subsidence amplitudes consistent with observed megathrust earthquakes. • Dynamic rupture simulations demonstrate physics-based controls on margin-wide rupture.