A document by Ryosuke Ando. Click on the document to view its contents.

Numerical simulations of Sequences of Earthquakes and Aseismic Slip (SEAS) have rapidly progressed to address fundamental problems in fault mechanics and provide self-consistent, physics-based frameworks to interpret and predict geophysical observations across spatial and temporal scales. To advance SEAS simulations with rigor and reproducibility, we pursue community efforts to verify numerical codes in an expanding suite of benchmarks. Here we present code comparison results from a new set of benchmark problems BP6-QD-A/S/C that consider a single aseismic slip transient induced by changes in pore fluid pressure consistent with fluid injection and diffusion in fault models with different treatments of fault friction. Ten modeling groups participated in problems BP6-QD-A and BP6-QD-S considering rate-and-state fault models using the aging and slip law formulations for frictional state evolution, respectively, allowing us to explore these ingredients across multiple codes and better understand how various computational factors affect the simulated evolution of pore pressure and aseismic slip. Comparisons of problems using the aging versus slip law illustrate how models of aseismic slip can differ in the timing and amount of slip achieved with different treatments of fault friction given the same perturbations in pore fluid pressure. We achieve excellent quantitative agreement across participating codes, with further agreement being found by ensuring sufficiently fine time-stepping and consistent treatment of remote boundary conditions. Our benchmark efforts offer a community-based example to reveal sensitivities of numerical modeling results, which is essential for advancing multi-physics SEAS models to better understand and construct reliable predictive models of fault dynamics.

Marine terraces have long been a subject of paleoseismology to reveal the rupture history of megathrust earthquakes. However, the mechanisms underlying their formation, in relation to crustal deformation, have not been adequately explained by kinematic models. A key challenge has been that the uplifted shoreline resulting from a megathrust earthquake tends to subside back to sea level during subsequent interseismic periods. This study focuses on the remaining permanent vertical deformation resulting from steady plate subduction and examines it quantitatively using three plate subduction models. Specifically, we pay attention to the effects of irregular geometries in the plate interface, such as subducted seamounts. Besides a simplified model examination, this study employs the plate geometry around the Sagami trough, central Japan, to compare with surface deformation observation. The subduction models employed are the kinematic subducting plate model, the elastic/viscoelastic fault model, and the mechanical subducting plate model (MSPM). The MSPM, introduced in this study, allows for more realistic simulations of crustal displacements by imposing net zero shear stress change on the plate boundary. Notably, the presence of a subducted seamount exerts a significant influence on surface deformation, resulting in a concentrated permanent uplift above it. The simulation of earthquake sequence demonstrates that coseismic uplifts can persist over time and contribute to the formation of marine terraces. The results demonstrated that the geological observations of coseismic and long-term deformations can be explained by the influence of a subducted seamount, previously identified in seismic surveys.

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 discovery of slow earthquakes illuminates the existence of a strange depth dependence of seismogenesis, which contradicts the common understanding of smooth brittle/seismic-ductile/aseismic transition as going deeper into the Earth’s surface layers. However, within the transitional layer on plate interfaces, observations have clarified slip velocities of slow earthquakes changing from those slower to faster with increasing depth, as described by the “seismogenic inversion layer.” We propose a new mechanical model that can consistently explain the classic brittle-ductile transition and this inversion phenomenon by considering the heterogeneous fault zone composed of brittle blocks in the ductile matrix. The key mechanism is the interplay between the volumetric fraction of brittle blocks and the viscosity of the surrounding plastically deformed matrix, where the former and the latter decrease with increasing temperature. This model is extended to shallow-slow earthquakes. Our results open a new pathway to infer the deformation mechanisms underlying slow earthquakes.

We developed a mechanical subducting plate model and re-examined the crustal deformation history in the Sagami Trough subduction zone, central Japan, the northernmost convergence boundary of the Philippine Sea Plate. The elevation distributions and formation ages of the Holocene marine terraces, representing past coseismic and long-term coastal uplifts, have been thoroughly investigated in this region. However, no physically consistent formation scenario to explain them has been demonstrated. Surface deformations within subduction zones are typically calculated using kinematic elastic dislocation models, such as the back-slip model. However, these models cannot explain permanent deformation after an earthquake sequence. This study develops a mechanical subducting plate model that balances the slips of interplate shear stress and can produce permanent deformations caused by a local bump geometry. We modeled earthquake recurrences by shear stress accumulation and coupling patches. As a result, we successfully reproduced the averaged uplift rate distribution estimated from the Holocene marine terraces. The findings suggest that the subducted seamount significantly affects long-term deformation patterns. In addition, the discrepancy between the elevation distributions and formation ages of Holocene marine terraces, which previous geological studies have indicated, can be interpreted by the rupture delay of coupling patches. This study also demonstrates that the traditional assumption of the back-slip model on the plate boundary for long-term subduction possibly results in an oversimplified model.

The 2008 Wenchuan Mw 7.9 mainshock has caused catastrophic destruction to cities along the northwestern margin of the Sichuan Basin. The Wenchuan-Maoxian Fault (WMF) on the hinterland side, along with a conjugate buried Lixian fault (LXF) was not activated by this earthquake but is likely to experience large earthquakes in the future. We perform 3D dynamic earthquake rupture simulations on the WMF and LXF to access the possibility of the earthquake occurrence and further explore the possible size of earthquakes and the distribution of high seismic risk in the future. We firstly invert focal mechanism solutions to get a heterogeneous tectonic stress field as the initial stress of simulation. Then we develop a new method to refine fault geometry through inverting long-term slip rates. Several fault nucleation points, friction coefficients, and initial stress states are tested, and the general rupture patterns for these earthquake scenarios are evaluated and could fall into three groups. Depending on initial conditions, the dynamic rupture may start in the LXF, leading to magnitude-7.0 earthquakes, or start in the WMF, then cascades through the LXF, leading to magnitude-7.5 earthquakes, or both start and arrest in the WMF, leading to around magnitude-6.5 or 7.0 earthquakes. We find that the rupture starting on the reverse oblique-slip tends to jump to the strike-slip fault, but the reverse process is suppressed. The rupture propagating eastward causes larger coseismic displacements than the westward propagation, and relatively high peak ground velocities are distributed near the northeastern end of WMF.