Seismic radiation and dynamic faulting have rarely been observed concurrently, limiting observational constraints on earthquake rupture processes and seismic hazards. Here we report the first directly validated observation of surface rupture and strong-motion generation captured by a security camera and corroborated by nearby seismic records and field measurements during the 2026 Kumamoto, Japan, earthquake. The near-fault ground motion comprised two velocity pulses generated by subsurface and surface fault ruptures, respectively.

Marine terraces have long been a subject of paleoseismology, revealing the rupture history of megathrust earthquakes. However, the crustal deformation mechanisms responsible for their formation remain inadequately explained by conventional kinematic models. A major challenge lies in the tendency of seismically uplifted shorelines to subside back to sea level during interseismic periods. This study focuses on the residual, permanent vertical deformation produced by repeated megathrust earthquakes. We investigate the effects of irregularities in the plate interface, particularly subducted seamounts. To address this, we introduce a mechanical subducting plate model (MSPM) that incorporates more realistic boundary conditions and three-dimensional geometry of the plate interface and subducting slab, using stress-boundary conditions. As a result, subducted seamounts significantly affect surface deformation, resulting in concentrated permanent uplift directly above them. We apply the MSPM to the geometry of the Sagami Trough, central Japan, and compare the simulation outcomes with the observations of marine terraces. The modeled earthquake sequences demonstrate that coseismic uplift can persist over time and contribute to terrace formation. These findings suggest that geological observations of both coseismic and long-term deformations can be explained by the influence of a subducted seamount, previously identified in seismic surveys.

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 quasi‐dynamic benchmark problems BP6‐QD‐A/S/C that consider an 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 (‐A) and slip (‐S) law formulations for frictional state evolution, respectively, allowing us to better understand how various computational factors across codes affect the simulated evolution of pore pressure and aseismic slip. Comparisons of problems using the aging versus slip law, and a constant friction coefficient (‐C), illustrate how aseismic slip models 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 attained by ensuring sufficiently fine time‐stepping and consistent treatment of 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.

This paper presents a fast H-matrices-based boundary integral equation method that can be applied to various wave-related problems. We propose an efficient algorithm for convolution in time direction for the intermediate domain between P- and S-waves using the exponentiation notation of the integral kernel. Furthermore, lattice H-matrices are used to address computational inefficiencies in parallel environments due to the hierarchical structure of H-matrices. 3D stress wave propagation simulations confirm the accuracy of the method.

The aging law and the slip law are two representative evolution laws of the rate- and state-dependent friction (RSF) law, based on canonical behaviours in three types of laboratory experiments: slide-hold-slide (SHS), velocity-step (VS) and steady-state (SS) tests. The aging law explains the SHS canon but contradicts the VS canon, and vice versa for the slip law. The later proposed composite law, which switches these two laws according to the slip rate V, explains both canons but contradicts the SS canon. This study constructs evolution laws satisfying all three canons throughout the range of variables where experiments have confirmed the canons. By recompiling the three canons, we have derived constraints on the evolution law and found that the evolution rates in the strengthening phases of the SHS and VS canons are so different that complete reconciliation throughout the entire range of variables is mathematically impossible. However, for the limited range of variables probed by experiments so far, we have found that the SHS and VS canons can be reconciled without violating the SS canon by switching the evolution function according to $\Omega$, the ratio of the state $\theta$ to its steady-state value $\theta _{\rm SS}$ for the instantaneous slip rate. We could generally show that, as long as the state evolution rate $\dot{\theta }$ depends only on the instantaneous values of V and $\theta$, simultaneous reproduction of the three canons, throughout the experimentally confirmed range, requires the aging-law-like evolution for $\Omega$ sufficiently below a threshold $\beta$ and the slip-law-like evolution for $\Omega$ sufficiently above $\beta$. The validity of the canons in existing experiments suggests $\beta \lesssim 0.01$.

The 2024 Mw 7.5 Noto Peninsula Earthquake broke through a previously documented active fault system over 150 km in the northern central Japanese Island. This fault system is characterized by geometrical complexity. It is important to understand the physical mechanism underlying the multi-fault rupture. We conduct fully dynamic rupture simulations and identify that the 3D fault geometry controls the observed rupture process and heterogeneous spatiotemporal patterns of the fault slip, seismic radiation and crustal deformation exhibiting about five meters of the maximum uplift. Aiming to examine the effect of the 3D fault geometry, we exclude the heterogeneity arising from the frictional properties. We also avoid retrospective frictional parameter tunings to fit the coseismic observations to test whether it is possible for our forward modeling to reproduce the coseismic observations. The 3D nonplanar geometry model is built based on the previously documented surface fault traces, and we use the regional stress field determined by the stress tensor inversion. As a result, the dynamic rupture simulation reasonably reproduces the observed characteristics of the heterogeneous deformation patterns. We find the rupture is accelerated, and slip is increased, where the fault is bent and optimally oriented to the regional stress orientations. Remarkably, the spatial distribution of surface displacement captured by the Synthetic Aperture Radar imageries is quantitatively reproduced, as characterized by two areas of large and small peaks of uplifts. Our findings may contribute to better constraining future earthquake rupture scenarios.

Earthquake-volcano interactions have been discussed to understand the underlying mechanisms of seismic ruptures or eruptions, yet the involvement of volcanic activity and the environment with fault slip termination remains unclear. Here, we present an unprecedented high-resolution image of fault motions and crustal structure at the rupture terminus in volcanic area from the 2016 Kumamoto earthquake by conducting synthetic aperture radar (SAR) data analysis and gravity inversion. We obtained a 3-D displacement field by applying multiple SAR analysis methods: standard SAR interferometry, split-bandwidth interferometry and pixel offset. We successfully mapped the ground displacements with a high-spatial resolution in the Aso caldera which was located on the eastern extension of the Futagawa fault that was the main source fault of this seismic event. We found that the rupture propagating on the Futagawa fault eastward penetrated into the Aso caldera and was divided into two major fault systems: a right-lateral fault system on the northern side and a left-lateral fault system on the southern side. However, they progressively converged immediately after penetrating into the caldera. A gravity-inferred 3-D density contrast structure revealed that a locally distributed low-density body existed in the shallow part (from the subsurface to a depth of similar to 3 km) of the western edge of the caldera. The slip distribution model showed that the slips on the bifurcated faults penetrated into the low-density region and subsequently dissipated. A numerical simulation on 3-D dynamic rupture demonstrated that the low-stress state in the caldera played a role in suppressing the rupture evolution. A thermally activated hydrothermal field has developed in the area where the fault slips were attenuated. We interpret that the hydrothermal system may create conditions favourable for low-stress field, and plastic properties in the hydrothermal environment may facilitate a further decrease in rock brittleness owing to the high temperature, resulting in the terminus of fault rupture.

The 2008 Wenchuan Mw 7.9 mainshock caused catastrophic destruction to cities along the northwestern margin of the Sichuan Basin. This earthquake did not activate the Wenchuan–Maoxian Fault (WMF) on the hinterland side and the conjugate buried Lixian Fault (LXF), but they could experience large earthquakes in the future. We propose a systematic scheme to develop scenario earthquakes for active fault systems with insufficient constrain of 3D fault geometries. We first performed stress tensor inversion to constrain the regional stress field. Then, we developed a new method to constrain fault geometries by inverting long-term slip rates under the given regional stress and applied it to the WMF. We conducted a set of 3D dynamic earthquake rupture simulations on the WMF and LXF to assess the scenarios of earthquake rupture processes. Several fault nucleation points, friction coefficients, and initial stress states are assessed, the general rupture patterns for these earthquake scenarios are evaluated, and finally, we find the scenarios that 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 cascade 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 jumps to the strike-slip fault, but the reverse process is impeded. Graphical Abstract

AbstractThe 2024 Mw 7.5 Noto Peninsula, Japan, earthquake was initiated within the source region of intense swarm activity. To reveal the mainshock early process, we relocated the earthquake hypocenters and found that many key phenomena, including the mainshock initiation, foreshocks, swarm earthquakes, and deep aseismic slip, occurred at parts of a previously unrecognized fault in intricate fault network. This fault is subparallel (several kilometers deeper) to a known active fault, and the mainshock initiation and foreshocks occurred at the front of a 2‐year westward swarm migration. The initiation location coincides with the destination of the upward migration of a deeper earthquake cluster via several smaller faults. Fluid supply, small earthquakes, and aseismic slip on the fault likely triggered the mainshock, leading to the first major rupture at the western region, propagating further to the west and east sides, resulting in an Mw7.5 event, exceeding 100 km in length.

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.

Large-scale earthquake sequence simulations using the boundary element method (BEM) incur extreme computational costs through multiplying a dense matrix with a slip rate vector. Hierarchical matrices (H-matrices) have often been used to accelerate this multiplication. However, the complexity of the structures of the H-matrices and the communication costs between processors limit their scalability, and they therefore cannot be used efficiently in distributed memory computer systems. Lattice H-matrices have recently been proposed as a tool to improve the parallel scalability of H-matrices. In this study, we developed a method for earthquake sequence simulations applicable to 3D nonplanar faults with lattice H-matrices. We present a simulation example and verify the mesh convergence of our method for a 3D nonplanar thrust fault using rectangular and triangular elements. We also performed performance and scalability analyses of our code. Our simulations, using over 10^5 degrees of freedom, demonstrated a parallel acceleration beyond 10^4 MPI processors and a >10-fold acceleration over the best performance when the normal H-matrices are used. Using this code, we can perform unprecedented large-scale earthquake sequence simulations on geometrically complex faults with supercomputers. The software HBI is made an open-source and freely available.

Numerical modeling of earthquake dynamics and derived insight for seismic hazard relies on credible, reproducible model results. The SEAS (Sequences of Earthquakes and Aseismic Slip) 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. Eight modeling groups participated in each benchmark, allowing us to explore these physical ingredients across multiple codes and better understand associated numerical considerations. We find that numerical resolution and computational domain size are critical parameters to obtain matching results, with increasing domain-size requirements posing challenges for volume-based codes even in 2D settings. Codes for BP1-FD implemented different criteria for switching between quasi-static and dynamic solvers, which require tuning to obtain matching results. In BP3-QD, proper remote boundaries conditions consistent with specified rigid body translation are required to obtain matching surface displacements. With these numerical and mathematical issues resolved, we obtain good agreement among codes in long-term fault behavior, earthquake recurrence intervals, and rupture features of peak slip rates and stress drops for both benchmarks. 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-sized characteristic earthquakes, and others exhibiting several earthquakes of differing magnitudes. These findings underscore the importance of considering full dynamics and non-vertical dip angles in SEAS models, as both influence short and long-term earthquake behavior, and associated hazards.

Fault bends are known to act as barriers to rupture propagation in many earthquakes. A recent compilation of the surface rupture traces of paleo earthquakes quantifies the probability of rupture termination as a function of the bend angle. To understand the physical basis of these unique statistics, we carry out 2D quasi-dynamic earthquake sequence simulations on a fault with either a restraining or releasing (double-) bend. The fault is loaded by steady sliding from an adjacent creeping region, rather than with backslip, together with a stress relaxation method to avoid unphysical stress buildup due to the curvature of faults. This relaxation is a proxy for unmodeled off-fault secondary faulting. This ensures the existence of a long-term steady earthquake cycle and stress field, the latter of which is determined by the balance between the slip-induced stress changes and relaxation terms. We quantify the influence of the bend on rupture propagation by computing passing ratio (i.e., the fraction of ruptures that propagate through the bend). Our simulations approximately reproduce the Biasi-Wesnousky empirical law for a wide range of parameters (e.g., bend width, stress relaxation time, background stress). Also, we reproduce geologic observations that the long-term slip rate of a fault has local minima at restraining bends, without assuming spatially varying loading rates using the backslip approach. Additionally, we find that restraining and releasing bends have different earthquake cycles. For restraining bends, many ruptures are arrested before reaching the center of the bend, and the passing ratio decreases with increasing the bend angle. For releasing bends, ruptures stop after passing the center of the bend, and the passing ratio abruptly drops from near unity to zero at around 20∘. This difference can be qualitatively explained by an energy balance approach. Our model has a potential to understand the seismogenesis of nonplanar faults and we can readily extend our model into 3D specific fault systems for hazard assessment.

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.

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.

We present a fast and memory-efficient algorithm for transient space-time-domain elastodynamic boundary-integral analysis. Associated data-sparse approximations and operations are named fast domain partitioning hierarchical matrices (FDP=H-matrices). The fast domain partitioning method (the FDPM) solves a known problem of hierarchical matrices (H-matrices) in compressing discretized elastodynamic kernel functions. A novel set of plane-wave approximations unites the FDPM and H-matrices in an accurate analytic manner. Memory usage is O(N log N) and computation time O(N M log N) in our algorithm for a single run with N boundary elements and M time steps. Consequent cost reduction is remarkable, considering the O((NM)-M-2) memory usage and O((NM2)-M-2) computational time to run the orthodox time-marching implementation. Numerical experiments verify FDP=H-matrices realize O(NM/log N) times smaller memory and computation time with ensuring the accuracy of integral analyses.

We quantitatively evaluated the emergence ages of the tectonically uplifted marine terraces, called the Numa terraces, in the southernmost part of the Boso Peninsula, central Japan. The combination of the complete dataset of the geological dating survey and the model inversion method is newly proposed. Along the Sagami Trough, the subduction boundary of the Philippine Sea Plate, M8 class interplate earthquakes are known to have repeatedly occurred and caused huge crustal deformation in adjacent regions, leaving tectonic landforms behind, such as the Numa terraces. Although many geological studies have investigated the Numa terraces, the estimation of their emergence ages contains ambiguity due to employed qualitative judgments and separate evaluations at different survey points. In this study, we first construct a unified dataset for dating throughout the southern Boso area by compiling the existing data and our new data. We also examin the lateral continuity of the terraces using a new technique to amplify subtle topographic changes in digital elevation data. Second, we propose a new method to estimate the emergence ages of terraces using the temporal distribution of the sample ages via a model of the sedimentation process and Bayesian statistic. The emergence ages of the three levels of the Numa terraces, uplifted in the prehistorical era, are determined as 5855–5455 yBP, 3345–3025 yBP, and 2125–1415 yBP, in descending order. This result robustly shows that the recurrence intervals of terrace-forming great Kanto earthquakes can vary by more than a factor of two. The robustness of the result is supported by the quantitative evaluations of the formation ages and the errors in the sedimentation process, made possible by our method for the first time. This result will provide inevitable information for further subduction zone researches and future hazard assessments.

In a dislocation problem, a paradoxical discordance is known to occur between an original smooth curve and an infinitesimally discretized curve. To solve this paradox, we have investigated a non-hypersingular expression for the integral kernel (called the stress Green’s function) which describes the stress field caused by the displacement discontinuity. We first develop a compact alternative expression of the non-hypersingular stress Green’s function for general 2-D and 3-D infinite homogeneous elastic media. We next compute the stress Green’s functions on a curved fault and revisit the paradox. We find that previously obtained non-hypersingular stress Green’s functions are incorrect for curved faults, and that smooth and infinitesimally segmented faults are equivalent. Their compatibility bridges the gap between analytical methods featuring curved faults and numerical methods using subdivided flat patches.