Seismic array processing is routinely used to infer detailed earthquake properties of intermediate and large events, however, the source properties of microseismicity often remain elusive. In this study, we use high signal-to-noise ratio seismograms of 204 earthquakes induced by the 6 km deep 2018 Espoo/Helsinki geothermal stimulation to evaluate the capabilities of beamforming and back-projection array methods. We show that mini array beamforming is sensitive to medium heterogeneity and requires calibration to mitigate systematic slowness biases.A combined and wave back-projection approach significantly improves depth resolution, reducing offsets to catalogue locations from km to m. Supported by numerical experiments, we demonstrate that back-projection swimming patterns can constrain focal mechanisms. Our results imply that back-projection of data collected over a wide azimuthal range can be used to monitor and characterize local-scale microseismicity, whereas beamforming calibration requires independently obtained reference observations.

We present 3-D spontaneous dynamic rupture earthquake scenarios for the Húsavík–Flatey Fault Zone (HFFZ) in Northern Iceland. We construct three fault system models consisting of up to 55 segments of varying geometric complexity. By varying hypocenter locations, we analyze rupture dynamics, fault interactions and their associated ground motions and observational uncertainties in 79 scenarios. We use regional observations to constrain 3-D subsurface velocities and viscoelastic attenuation as well as fault stress and strength. Our models account for topo-bathymetry, off-fault plasticity and we explore the effect of fault roughness. Our spontaneous dynamic rupture scenarios can match historic magnitudes. We show that the fault system segmentation and geometry, hypocenter locations, initial stress conditions and fault roughness have strong effects on multi-fault rupture dynamics across the HFFZ. Breaking of different portions of the same fault system leads to varying rupture dynamics, slip distributions and magnitudes. All dynamic rupture scenarios yield highly heterogeneous near-field ground motions. We observe amplification from rupture directivity, geometric complexities, and amplification and shielding due to topography. We recover a magnitude-consistent attenuation relationship in good agreement with new regional empirical ground motion models. Physics-based ground motion variability changes with distance and increases for unilateral vs. bilateral rupture. Our study illustrates important ingredients for fully physics-based, regional earthquake scenarios, their respective importance for rupture dynamics and ground motion modeling and how they can be observationally constrained and verified. We entail that dynamic rupture scenarios can be useful for non-ergodic probabilistic seismic hazard assessment, specifically in data-limited regions.

We study the effects of pore fluid pressure (Pf ) on the pre-earthquake, near-fault stress state and 3D earthquake rupture dynamics through 6 scenarios utilizing a structural model based on the 2004 Mw 9.1 Sumatra-Andaman earthquake. As pre-earthquake Pf magnitude increases, effective normal stress and fault shear strength decrease. As a result, magnitude, slip, peak slip rate, stress drop and rupture velocity of the scenario earthquakes decrease. Comparison of results with observations of the 2004 earthquake support that pre-earthquake Pf averages near 97 % of lithostatic pressure, leading to pre-earthquake average shear and effective normal tractions of 4-5 MPa and 22 MPa. The megathrust in these scenarios is weak, in terms of low mean shear traction at static failure and low dynamic friction coefficient during rupture. Apparent co-seismic principal stress rotations and absolute post-seismic stresses in these scenarios are consistent with the variety of observed aftershock focal mechanisms. In all scenarios, the mean apparent stress rotations are larger above than below the megathrust. Scenarios with larger Pf magnitudes exhibit lower mean apparent principal stress rotations. We further evaluate pre-earthquake Pf depth distribution. If Pf follows a sublithostatic gradient, pre-earthquake effective normal stress increases with depth. If Pf follows the lithostatic gradient exactly, then this normal stress is constant, shifting peak slip and peak slip rate up-dip. This renders constraints on near-trench strength and constitutive behavior crucial for mitigating hazard. These scenarios provide opportunity for future calibration with site-specific measurements to constrain dynamically plausible megathrust strength and Pf gradients.

From interpreting data to scenario modeling of subduction events, numerical modeling has been crucial for studying tsunami generation by earthquakes. Seafloor instruments in the source region feature complex signals containing a superposition of seismic, ocean acoustic, and tsunami waves. Rigorous modeling is required to interpret these data and use them for tsunami early warning. However, previous studies utilize separate earthquake and tsunami models, with one-way coupling between them and approximations that might limit the applicability of the modeling technique. In this study, we compare four earthquake-tsunami modeling techniques, highlighting assumptions that affect the results, and discuss which techniques are appropriate for various applications. Most techniques couple a 3D Earth model with a 2D depth-averaged shallow water tsunami model. Assuming the ocean is incompressible and that tsunami propagation is negligible over the earthquake duration leads to technique (1), which equates earthquake seafloor uplift to initial tsunami sea surface height. For longer duration earthquakes, it is appropriate to follow technique (2), which uses time-dependent earthquake seafloor velocity as a time-dependent forcing in the tsunami mass balance. Neither technique captures ocean acoustic waves, motivating newer techniques that capture the seismic and ocean acoustic response as well as tsunamis. Saito et al. (2019) propose technique (3), which solves the 3D elastic and acoustic equations to model the earthquake rupture, seismic wavefield, and response of a compressible ocean without gravity. Then, sea surface height is used as a forcing term in a tsunami simulation. A superposition of the earthquake and tsunami solutions provides the complete wavefield, with one-way coupling. The complete wavefield is also captured in technique (4), which utilizes a fully-coupled solid Earth and ocean model with gravity (Lotto & Dunham, 2015). This technique, recently incorporated into the 3D code SeisSol, simultaneously solves earthquake rupture, seismic waves, and ocean response (including gravity). Furthermore, we show how technique (3) follows from (4) subject to well-justified approximations.

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.

Despite decades-long debate over the mechanics of low-angle normal faults dipping less than 30°, many questions about their strength, stress, and slip remain unresolved. Recent geologic and geophysical observations have confirmed that gently-dipping detachment faults can slip at such shallow dips and host moderate-to-large earthquakes. Here, we analyze the first 3D dynamic rupture models to assess how different stress and strength conditions affect rupture characteristics of low-angle normal fault earthquakes. We model observationally constrained spontaneous rupture under different loading conditions on the active Mai’iu fault in Papua New Guinea, which dips 16-24° at the surface and accommodates ~8 mm/yr of horizontal extension. We analyze four distinct fault-local stress scenarios: 1) Andersonian extension, as inferred in the hanging wall; 2) back-rotated principal stresses inferred paleopiezometrically from the exhumed footwall; 3) favorably rotated principal stresses well-aligned for low-angle normal-sense slip; and 4) Andersonian extension derived from depth-variable static fault friction decreasing towards the surface. Our modeling suggests that subcritically stressed detachment faults can host moderate earthquakes within purely Andersonian stress fields. Near-surface rupture is impeded by free-surface stress interactions and dynamic effects of the gently-dipping geometry and frictionally stable gouges of the shallowest portion of the fault. Although favorably-inclined principal stresses have been proposed for some detachments, these conditions are not necessary for seismic slip on these faults. Our results demonstrate how integrated geophysical and geologic observations can constrain dynamic rupture model parameters to develop realistic rupture scenarios of active faults that may pose significant seismic and tsunami hazards to nearby communities.

The theory of rate and state friction unifies field, laboratory, and theoretical analysis of the evolution of slip on natural faults. While the observational study of earthquakes and aseismic fault slip is hampered by its strong multi-scale character in space and time, numerical simulations are well-positioned to link the laboratory study of grain-scale processes to the scale at which rock masses move. However, challenges remain in accurately representing the complex and permanently evolving sub-surface fault networks that exist in nature. Additionally, the common representation of faults as interfaces may miss important physical aspects governing volumetric fault system behavior. In response, we propose a transient viscous rheology that produces shear bands that closely mimic the rate- and state-dependent sliding behavior of equivalent fault interfaces. Critically, we show that the expected tendency of the continuum rheology for runaway localization and mesh dependence can be halted by including an artificial diffusion-type regularization of anelastic strain rate in the softening law. We demonstrate analytically and numerically using a simplified fault transect that important aspects of the frictional behavior are not significantly affected by the introduced regularization. Any discrepancies with respect to the interfacial description of fault behavior are critically evaluated using 1D numerical velocity stepping and spring-slider experiments. ;Since no new physical parameters are introduced, our model may be straightforwardly used to extend the existing modeling techniques. The model predicts the emergence of complex patterns of shear localization and delocalization that may inform the interpretation of complex damage distributions observed around faults in nature.

Advances in physics-based earthquake simulations, utilizing high-performance computing, have been exploited to better understand the generation and characteristics of the high-frequency seismic wavefield. However, direct comparison to ground motion observations of a specific earthquake is challenging. We here propose a new approach to simulate data-fused broadband ground motion synthetics using 3D dynamic rupture modeling of the 2016 Mw6.2 Amatrice, Italy earthquake. We augment a smooth, best-fitting model from Bayesian dynamic rupture source inversion of strong-motion data (<1 Hz) with fractal fault roughness, frictional heterogeneities, viscoelastic attenuation, and topography. The required consistency at long periods allows us to quantify the role of dynamic source heterogeneities, such as the 3D roughness drag, from observational broadband seismic waveforms. We demonstrate that 3D data-constrained fully dynamic rupture synthetics show good agreement with various observed ground-motion metrics up to ~5 Hz and are an important avenue towards non-ergodic, physics-based seismic hazard assessment.