Dynamic restrengthening and fault heterogeneity explain megathrust earthquake complexity

Abstract

Megathrusts host Earth's largest earthquakes. Understanding the physical conditions controlling their rupture dynamics is critical for assessing seismic and tsunami hazards. These earthquakes often display complex rupture dynamics, exemplified by the 2011 Tohoku-Oki earthquake, which exhibited multiple rupture episodes, depth-dependent seismic radiation, and substantial tsunamigenic slip near the trench. However, how such complexity arises from preexisting physical conditions remains uncertain. Here, we demonstrate that the observed rupture complexity of the Tohoku-Oki earthquake can spontaneously and self-consistently emerge, driven by rapid coseismic frictional restrengthening and data-informed fault heterogeneity. We use an ensemble of 3D dynamic rupture simulations to identify that mixed downdip pulse-like and updip crack-like rupture are driven by dynamic stress redistribution with episodic rupture reactivation. By featuring low fault strength compared to its dynamic stress drop, a preferred model can consistently reproduce the observed complex depth-dependent propagation speeds, multiple rupture fronts as imaged by back-projection, and large tsunamigenic slip at the trench. Our findings demonstrate that preexisting fault heterogeneity conjointly with dynamic frictional weakening and restrengthening drives seemingly unexpected megathrust rupture complexity, highlighting the need to include dynamic effects into physics-based seismic and tsunami hazard assessments of future earthquakes.

Figures (39)

• Figure 1: Overview of 3D dynamic rupture simulations of the 2011 Tohoku-Oki earthquake and their initial conditions. Depth contours (gray, 10 km intervals) and hypocenter location (star) are shown in all panels. (a) Snapshot of the simulated absolute slip rate and seismic wavefield evolution (vertical particle velocity) at 66 s, highlighting multiple reactivated slip pulses propagating downdip and crack-like rupture accumulating large slip near the trench. The model incorporates realistic slab geometry and high-resolution topobathymetry within an unstructured tetrahedral mesh refined near the slab and onshore region. (b) Depth-dependent frictional properties ($a-b$, see Methods Sec. "\ref{['Fault friction']}") with velocity-strengthening behavior in the shallow ($<$9 km) and deep ($>$45 km) regions, transitioning to velocity-weakening in the seismogenic zone. (c) Laterally homogeneous, depth-dependent ambient stress and frictional strength initial conditions informed by regional tectonics (see Methods Sec. "\ref{['Prestress']}"), showing the relative prestress ratio $R$ (maximum possible stress drop over frictional strength drop, Eq.\ref{['eq:R']}, Aochi20031999). The principal stress direction ($\sigma_1$ at an azimuth of 100$^\circ$ and a plunge angle of 8$^\circ$, Heidbach2018World) is indicated with arrows. (d) Heterogeneous stress initial conditions combining the ambient background stress shown in (c) and heterogeneous initial stress inferred from the median slip distribution of 32 data-constrained slip models (Supplementary Figs. \ref{['SFig:Kinematic_med_slip_model']}, \ref{['SFig:Traction_cij']}, Wong2024Quantitative). Along-dip initial shear stresses are shown in Supplementary Fig. \ref{['FigE1:Initial stress conditions']}.
• Figure 2: Preferred 3D dynamic rupture scenario of the Tohoku-Oki earthquake constrained by geodetic observations and seismic moment release rate. Gray contours indicate depth (10 km intervals), and the star is the hypocenter location Hayes2011Rapid. Rupture extends 200 km along-dip and 360 km along-strike, producing a moment magnitude of $M_W~8.97$ and a duration of 120 s. The total radiated seismic energy is $\approx 7.7\times 10^{17}J$, within observational estimates of $4.2-9.1 \times 10^{17}J$ for the Tohoku-Oki earthquake Ide2011ShallowLay2012Depthvarying. (a) Fault slip distribution with comparison between observed and simulated geodetic displacements onshore and offshore. Black arrows denote observed horizontal displacements from offshore and onshore stations. Blue and red arrows represent simulated horizontal displacements offshore and onshore, respectively, achieving variance reductions of 77% (onshore) and 55% (offshore). (b) Synthetic moment rate release compared with observational inferences from teleseismic by the USGS model Hayes2011Rapid and SCARDEC inversion results VallA2016New. Heterogeneous spatial distributions of (c) peak slip rate, (d) rupture speed, (e) rupture front timing (10 s intervals, gray contours), and (f) along-dip stress drop.
• Figure 3: Repeated dynamic rupture reactivation enabled by rapid coseismic weakening and restrengthening during the preferred Tohoku-Oki earthquake dynamic rupture model. (a) Map-view snapshots of rupture evolution from 10 s to 85 s simulation time, showing three main re-nucleation episodes at 15 s, between 25--40 s, and at 50 s. The white star indicates the hypocenter. Similarly "spiraling" rupture fronts have been observed in recent laboratory experiments Cochard2024Propagation. Supplementary Fig. \ref{['FigE2:All_slip_rate_5s']} and Supplementary Video S1 show the complete rupture evolution. Supplementary Fig. \ref{['FigE3:Spiral_rupture']} shows the detailed evolution of "spiraling" rupture fronts. (b) Slip rate evolution along a dip profile through the hypocenter, highlighting multiple episodes of rupture reactivation. Crosses and circles indicate the locations of high-frequency radiation from back-projection using the US and European arrays, respectively Meng2011Window. Rupture propagates faster updip ($\approx$2.5 km/s) compared to downdip ($\approx$1.7 km/s), matching the observational results from the back-projection analyses Meng2011Window. (c)-(d) Temporal evolution of along-dip shear stress (purple) and effective friction coefficient (blue, Methods " \ref{['Fault friction']}") along the hypocentral dip profile, highlighting rapid variations coincident with dynamic rupture reactivation; lighter colors indicate higher values. (e) Time series at the hypocenter of slip rate (red), along-dip shear stress (purple), and effective friction coefficient (blue), showing repeated rupture reactivation (slip rate $\ge$ 0.05 m/s, shaded red) and rapid frictional restrengthening. Supplementary Fig. \ref{['FigE3:Profile_comp']} shows updip and downdip time series evolution of slip rate, along-dip shear stress, and effective friction coefficient.
• Figure 4: Depth-varying rupture styles featuring downdip short-duration slip pulses and updip large-slip crack-like ruptures. (a) Slip-rate evolution with depth along the same hypocenteral dip profile in Fig. \ref{['Fig3:reactivation']}, highlighting crack-like ruptures shallower than 15 km depth (red) and pulse-like ruptures at greater depths (blue). (b) Power spectra of the slip rates, illustrating systematically higher-frequency content in downdip pulses compared to the shallower crack-like ruptures. The shallow crack-like rupture spectra follow a -1 slope over the 0.02-1 Hz range, whereas the pulse-like rupture exhibits a shallower -0.7 slope within the same frequency band. (c--d) Spatial distribution of high-pass filtered (HP) and low-pass filtered (LP) peak slip rates, respectively. (e) Spatial variation in the ratio of high-frequency to low-frequency peak slip rates (HP/LP). Shallow regions ($<$15 km depth) exhibit predominantly crack-like rupture with low HP/LP ratios, while downdip and hypocentral areas show pulse-like ruptures enriched in high-frequency content, consistent with observations from back-projection and regional strong-ground motion analyses Meng2011WindowKurahashi2013ShortPeriod. Supplementary Fig. \ref{['FigE4:Depth_dependent_freq_ranges']} shows additional comparisons in three frequency ranges: 10--2 s, 2--0.5 s, and $>$ 0.5 s.
• Figure 5: Rapid transverse expansion of a circular rupture front, resembling recent laboratory observations Cochard2024Propagation. Snapshots show slip rate evolution of a circular, spiraling rupture front from 64 s to 69 s. The black arrows mark the onset location where the circular rupture front forms. The spiraling rupture reaches a radius of approximately 12 km from the hypocenter and sweeps a 180$^{\circ}$ arc within 5 s. Its transverse rupture front expansion speed exceeds 7 km/s, which is higher than the local shear wave speed at the corresponding depth of $\approx$20 km (Supplementary Table S1).