The September 2018, M w 7.5 Sulawesi earthquake occurring on the Palu-Koro strike-slip fault system was followed by an unexpected localized tsunami. We show that direct earthquakeinduced uplift and subsidence could have sourced the observed tsunami within Palu Bay. To this end, we use a physics-based, coupled earthquake-tsunami modeling framework tightly constrained by observations. The model combines rupture dynamics, seismic wave propagation, tsunami propagation and inundation. The earthquake scenario, featuring sustained supershear rupture propagation, matches key observed earthquake characteristics, including the moment magnitude, rupture duration, fault plane solution, teleseismic waveforms and inferred horizontal ground displacements. The remote stress regime reflecting regional transtension applied in the model produces a combination of up to 6 m left-lateral slip and up to 2 m normal slip on the straight fault segment dipping 65 East beneath Palu Bay. The time-dependent, 3D seafloor displacements are translated into bathymetry perturbations with a mean vertical offset of 1.5 m across the submarine fault segment. This sources a tsunami with wave amplitudes and periods that match those measured at the Pantoloan wave gauge and inundation that reproduces observations from field surveys. We conclude that a source related to earthquake displacements is probable and that landsliding may not have been the primary source of the tsunami. These results have important implications for submarine strike-slip fault systems worldwide. Physics-based modeling offers rapid response specifically in tectonic settings that are currently underrepresented in operational tsunami hazard assessment.

We present a high-resolution simulation of the 2004 Sumatra-Andaman earthquake, including non-linear frictional failure on a megathrustsplay fault system. Our method exploits unstructured meshes capturing the complicated geometries in subduction zones that are crucial to understand large earthquakes and tsunami generation. These up-to-date largest and longest dynamic rupture simulations enable analysis of dynamic source effects on the seafloor displacements.To tackle the extreme size of this scenario an end-to-end optimization of the simulation code SeisSol was necessary. We implemented a new cache-aware wave propagation scheme and optimized the dynamic rupture kernels using code generation. We established a novel clustered local-time-stepping scheme for dynamic rupture. In total, we achieved a speed-up of 13.6 compared to the previous implementation. For the Sumatra scenario with 221 million elements this reduced the time-to-solution to 13.9 hours on 86,016 Haswell cores. Furthermore, we used asynchronous output to overlap I/O and compute time.

Contractional sandbox experiments that simulate crustal accretion and direct shear tests both provide direct data on the amount of work required to create faults (W prop) in granular materials. Measurements of force changes associated with faulting reveal the work consumed by fault growth, which can be used to predict fault growth path and timing. Within the contractional experiments, the sequence and style of early faulting is consistent for the range of sand pack thicknesses tested, from 12 to 30 mm. Contrary to expectations that W prop is only a material property, the experimental data show that for the same material, W prop increases with sand pack thickness. This normal stress dependence stems from the frictional nature of granular materials. With the same static and sliding friction values, incipient faults initiated deeper in the sand pack have larger shear stress drops, due to increased normal compression, σ n. For CV32 sand, the relationship between W prop and σ n, calculated from the force drop data as W prop (J/m 2) = 2.0x10-4 (m) σ n (Pa), is consistent with the relationship calculated from direct shear test data as W prop (J/m 2) = 2.4x10-4 (m) σ n (Pa). Testing of different materials within the contractional sandbox (fine sand and glass beads) shows the sensitivity of W prop to material properties. Both material properties and normal stress should be considered in calculations of the work consumed by fault growth in both analog experiments and crustal fault systems.

No abstract

Taking the full complexity of subduction zones into account is important for realistic modeling and hazard assessment of subduction zone seismicity and associated tsunamis. Studying seismicity requires numerical methods that span a large range of spatial and temporal scales. We present the first coupled framework that resolves subduction dynamics over millions of years and earthquake dynamics down to fractions of a second. Using a two-dimensional geodynamic seismic cycle (SC) model, we model 4 million years of subduction followed by cycles of spontaneous megathrust events. At the initiation of one such SC event, we export the self-consistent fault and surface geometry, fault stress and strength, and heterogeneous material properties to a dynamic rupture (DR) model. Coupling leads to spontaneous dynamic rupture nucleation, propagation, and arrest with the same spatial characteristics as in the SC model. It also results in a similar material-dependent stress drop, although dynamic slip is significantly larger. The DR event shows a high degree of complexity, featuring various rupture styles and speeds, precursory phases, and fault reactivation. Compared to a coupled model with homogeneous material properties, accounting for realistic lithological contrasts doubles the amount of maximum slip, introduces local pulse-like rupture episodes, and relocates the peak slip from near the downdip limit of the seismogenic zone to the updip limit. When an SC splay fault is included in the DR model, the rupture prefers the splay over the shallow megathrust, although wave reflections do activate the megathrust afterward.

Megathrust faults host the largest earthquakes on Earth which can trigger cascading hazards such as devastating tsunamis. Determining characteristics that control subduction zone earthquake and tsunami dynamics is critical to mitigate megathrust hazards, but is impeded by structural complexity, large spatio-temporal scales, and scarce or asymmetric instrumental coverage. Here we use high-performance computing multi-physics simulations to show that tsunamigenesis and earthquake dynamics are controlled by along-arc variability in regional tectonic stresses together with depth-dependent variations in rigidity and yield strength of near-fault sediments. We aim to identify dominant regional factors controlling megathrust hazards. To this end, we demonstrate how to unify and verify the required initial conditions for geometrically complex, multi-physics earthquake-tsunami modeling from interdisciplinary geophysical observations. We present large-scale computational models of the 2004 Sumatra-Andaman earthquake and Indian Ocean tsunami that reconcile near-and far-field seismic, geodetic, geological, and tsunami observations and reveal tsunamigenic trade-offs between slip to the trench, splay faulting, and bulk yielding of the accretionary wedge. Our computational capabilities render possible the incorporation of present and emerging high-resolution observations into dynamic-rupture-tsunami models and will be applicable to other large megathrust earthquakes. Our findings highlight the importance of regional-scale structural heterogeneity to decipher megathrust hazards.Variations in megathrust earthquake rupture behaviour are associated with tectonic, mechanical and structural factors highlighting the importance of depth-dependent and along-arc subduction zone heterogeneity (1-5). Large tsunamis may be caused by various co-seismic mechanisms including large slip to the trench, as observed during the 2011 Tohoku earthquake (6) and inferred from

Taking the full complexity of subduction zones into account is important for realistic modelling and hazard assessment of subduction zone seismicity and associated tsunamis. Studying seismicity requires numerical methods that span a large range of spatial and temporal scales. We present the first coupled framework that resolves subduction dynamics over millions of years and earthquake dynamics down to fractions of a second. Using a two-dimensional geodynamic seismic cycle (SC) model, we model 4~million years of subduction followed by cycles of spontaneous megathrust events. At the initiation of one such SC event, we export the self-consistent fault and surface geometry, fault stress and strength, and heterogeneous material properties to a dynamic rupture (DR) model. Coupling leads to spontaneous dynamic rupture nucleation, propagation and arrest with the same spatial characteristics as in the SC model. It also results in a similar material-dependent stress drop, although dynamic slip is significantly larger. The DR event shows a high degree of complexity, featuring various rupture styles and speeds, precursory phases, and fault reactivation. Compared to a coupled model with homogeneous material properties, accounting for realistic lithological contrasts doubles the amount of maximum slip, introduces local pulse-like rupture episodes, and relocates the peak slip from near the downdip limit of the seismogenic zone to the updip limit. When an SC splay fault is included in the DR model, the rupture prefers the splay over the shallow megathrust, although wave reflections do activate the megathrust afterwards.

Restraining bends along strike-slip fault systems evolve by both propagation of new faults and abandonment of fault segments. Scaled analog modeling using wet kaolin allows for qualitative and quantitative observations of this evolution. To explore how bend geometry affects evolution, we model bends with a variety of initial angles, θ, from θ = 0°for a straight fault to θ = 30°. High-angle restraining bends (θ ≥ 20°) overcome initial inefficiencies by abandoning unfavorably oriented restraining segments and propagating multiple new, inwardly dipping, oblique-slip faults that are well oriented to accommodate convergence within the bend. Restraining bends with 0°< θ ≤ 15°maintain activity along the restraining bend segment and grow a single new oblique slip fault on one side of the bend. In all restraining bends, the first new fault propagates at~5 mm of accumulated convergence. Particle Image Velocimetry analysis provides a complete velocity field throughout the experiments. From these data, we quantify the strike-slip efficiency of the system as the percentage of applied plate-parallel velocity accommodated as slip in the direction of plate motion along faults within the restraining bend. Bends with small θ initially have higher strike-slip efficiency compared to bends with large θ. Although they have different fault geometries, all systems with a 5 cm bend width reach a steady strike-slip efficiency of 80% after 50 mm of applied plate displacement. These experimental restraining bends resemble crustal faults in their asymmetric fault growth, asymmetric topographic gradient, and strike-slip efficiency.