Search for: All records

Abstract Understanding earthquake dynamics is essential for seismic hazard assessment and risk mitigation. In this context, Bayesian inference provides valuable insights into model parameters by combining simulation models with real-world data. Such Bayesian parameter inference with uncertainty quantification (UQ) requires numerous simulation runs and is therefore often computationally infeasible. Already, a single high-fidelity earthquake simulation – governed by a linear hyperbolic seismic wave equation coupled nonlinearly to a friction law and plastic deformation – is computationally expensive. In this study, we investigate the use of fused ensemble simulations to accelerate large earthquake simulation workflows and UQ studies. We implement and evaluate this approach in SeisSol, a high-performance computing software for the simulation of complex earthquake events that uses an Arbitrary high-order DERivative Discontinuous Galerkin (ADER-DG) scheme. Via fused ensembles, we turn the element-local small sparse/dense matrix operations into tensor contractions working on a dense rank-3 tensor and sparse matrices. These are again executed via loops of small, sparse/dense matrix operations, but with better computational efficiency, due to better exploitation of SIMD instructions on CPUs. We also compare two implementation variants (with different implementation effort) for kernels modeling non-linear dynamic rupture and effects of material plasticity. Our results demonstrate that fused simulations can be up to 5.54 times faster than a single execution – though this depends strongly on the discretization order, the problem size, and the compute architecture. For a full UQ example workflow, we demonstrate a speedup of 1.6, resulting in 36% savings in node hours for the entire workflow. 

Abstract Slow slip events (SSEs) in the Cascadia subduction zone exhibit along‐strike segmentation, where the central segment has longer recurrence intervals but smaller moments. We quantify the controls on this variability by combining geodetic inter‐SSE coupling inversion with Bayesian inference of a quasi‐dynamic rate‐and‐state friction SSE‐cycle model accelerated by reduced‐order modeling. Our simulations show that effective normal stress controls the SSE recurrence interval, while subduction coupling controls the SSE moment. Our inversion of inter‐SSE GNSS velocities yields comparable inter‐SSE coupling along strike (∼60%–70%), whereas long‐term coupling ranges from near zero in the south to ∼42% in central Cascadia. Transient SSEs recover the full slip deficit in the south but leave persistent deficits of ∼30% and ∼42% of plate convergence in the north and central segments. These results indicate that effective normal stress and inter‐SSE coupling provide a unified geodetic and physics‐based explanation for Cascadia SSE segmentation. 

Abstract Induced seismicity associated with subsurface fluid injection poses significant challenges for geothermal reservoir management. While quasi-static models capture long-term pressure, temperature, and stress evolution, dynamic rupture codes resolve transient earthquake mechanics. In this study, we develop a fully coupled framework that integrates both approaches to simulate fault reactivation and induced seismicity in a synthetic geothermal system. Unlike previous one-way coupled frameworks, our approach implements a closed-loop feedback, where stress redistribution from dynamically simulated seismic events is fed back into the quasi-static THM reservoir model. A 3D thermo-hydro-mechanical model tracks stress build-up during injection until a failure criterion is met, triggering a transition to a dynamic rupture solver that computes spontaneous 3D rupture propagation and seismic wave radiation. An automated coupling protocol ensures seamless data exchange between solvers. We demonstrate the framework by simulating two injection-induced events on a single-fault reservoir model. Results show that event-driven fault displacements and dynamic stress transfers significantly alter reservoir stress states. Simulated Coulomb stress changes, focal mechanisms, and slip distributions align with the ambient stress field and fault geometry. This closed-loop approach provides a robust tool for assessing seismic risk and optimizing injection strategies in geothermal environments. 

ABSTRACT Recent advances in Probabilistic Seismic Hazard Analysis (PSHA) leverage physics-based ground-motion simulations to estimate seismic hazard, such as the CyberShake project. However, computational costs quickly escalate when performing PSHA for numerous faults or sites and can become prohibitively expensive. To reduce computational demands, CyberShake uses reciprocity and interpolates physics-informed corrections from simulations conducted at fewer locations, but the accuracy of these interpolations remains poorly quantified. To quantify the interpolation accuracy, we derive high-resolution, frequency-dependent site terms for southern California and compare them with interpolated site terms using the CyberShake approach. We accomplish this by performing a set of earthquake point-source simulations distributed along the nonplanar fault geometry for the southern San Andreas fault (SSAF) extending from Bombay Beach to Lake Hughes. Using SeisSol, we simulate three minutes of viscoelastic seismic wave propagation for these sources and store the horizontal-component Green’s functions for 480,000 sites. We then use a scientific machine learning approach based on interpolated proper orthogonal decomposition to construct an accurate reduced-order model of the Green’s functions to efficiently predict effective amplitude spectra (EAS) for finite-source rupture models of SSAF earthquakes. Using minimum curvature interpolation with tension, as used in CyberShake, we compare the interpolated site terms against our high-resolution site terms. We identify local discrepancies with EAS differing by up to a factor of approximately three. Furthermore, we identify locations where unexpectedly high or low ground motions are missed when using the interpolated dataset for these earthquakes. We estimate that our approach may be used within CyberShake to reduce the time-to-solution by a factor of 336 for the entire earthquake rupture forecast. Our analysis of physics-based site terms provides more insight into the seismic hazard due to SSAF ruptures and guides future developments by combining high-performance computing and reduced-order modeling techniques for PSHA. 

Abstract Simulating sequences of earthquakes and aseismic slip (SEAS) on realistic fault systems provides a physical framework to understand the evolution of the seismic cycle, but remains computationally expensive. Volumetric approaches offer the physical flexibility to handle complex geometries and heterogeneous off-fault media but may incur prohibitively high computational costs when applied to the vast range of spatial and temporal scales inherent to earthquake cycles. This article documents developments in Tandem, an open-source volumetric SEAS simulation software package that addresses these challenges using a symmetric interior penalty discontinuous Galerkin formulation on unstructured curvilinear meshes in 2D and 3D with high-order polynomial basis functions. We describe Tandem from a user’s perspective, covering mesh formats, checkpointing, model configuration via human-readable files, and flexible loading schemes for various tectonic settings. To lower barriers to usage, Tandem is distributed as a standalone C++/PETSc code, a preconfigured virtual machine image, and as an application on the Quakeworx Science Gateway, where users can run simulations without installing dependencies or securing direct high-performance computing (HPC) access. We report on practical development choices of interest to SEAS and other scientific software developers. To mitigate the cost of volumetric discretization, the software provides both matrix-free and assembled-matrix formulations with a fully volumetric explicit approach utilizing hybrid geometric–algebraic multigrid preconditioners and a discrete Green’s function mode that accelerates timestepping by precomputing traction kernels. We document performance on both central processing unit and graphical processing unit supercomputers, including near-ideal weak scaling up to 112,000 message passing interface (MPI) ranks on the supercomputer Frontera. Finally, as demonstration examples, we present validation in a 3D SEAS community benchmark (BP7) and a 2D example examining how off-megathrust material heterogeneity influences seismic cycle behavior. By sharing software design choices as well as practical guidance for its use, we hope to make volumetric HPC-driven SEAS modeling more accessible to the earthquake science community. 

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. 

Abstract Understanding the mechanisms controlling deformation localization is crucial for our understanding of fault mechanics and improving seismic hazard assessment, but has not been extensively studied for normal‐faulting earthquakes. Here, we present a thorough analysis of the 2016 Mw 6.5 Norcia, Italy, earthquake using high‐resolution satellite geodesy. We investigate the degree of deformation localization, evaluate its controlling factors and the link with the distribution of slip with depth. Using the optical image correlation technique with an innovative method for noise correction, we measure the near‐field 3D displacements associated with the Norcia event. Based on these measurements, we quantify the amount of off‐fault deformation (OFD) and evaluate how it varies with external factors. We find 46 (25 cm) of OFD on average and a strong correlation with local topographic slope and near‐surface lithology, with increased distributed deformation where ruptures traverse unconsolidated sediments and areas of gentler slopes. In contrast, the correlation between OFD and the fault segment orientation relative to the regional stress field is weak. We develop a comprehensive slip model accounting for complex multi‐segmented fault geometry, topography, and 3D elastic structure through a joint inversion of optical, InSAR, and GPS data. The inversion reveals a highly heterogeneous slip distribution characterized by large slip (up to 3.5 m) at depth, and several shallow slip patches. We find a pronounced average shallow slip deficit (SSD) of 73, with no along‐strike correlation between SSD and OFD. This suggests that the OFD primarily reflects surficial inelastic processes occurring within the shallow soil. 

SUMMARY The Hellenic Arc Subduction Zone is the most seismically active region in the Mediterranean, capable of generating large earthquakes and tsunami. Given the proximity of densely populated coastlines, understanding the characteristics of potential future large tsunamigenic earthquakes is crucial for assessing tsunami hazard. We present nonlinear shallow water tsunami simulations sourced from the static seafloor displacements of five $$M_{\rm W}$$ 8–9 dynamic rupture earthquake scenarios along the Hellenic Arc, varying in hypocentral location, rupture extent and moment magnitude. Our results show that shallow slip and off-fault plastic deformation control the vertical near-trench uplift, tsunami height, and tsunami arrival time, particularly towards the south-west of the Hellenic Arc. Plasticity increases the tsunami amplitudes by over 50 percent. Maximum tsunami amplitudes reach up to 6.6 m near central and eastern Mediterranean coastlines, while the northern Aegean and the western Mediterranean remain mostly shielded by landmasses in all scenarios. We use the time-dependent seafloor deformation of the worst-case scenario among those considered to highlight differences between static and time-dependent sourcing for tsunami modelling. The results show up to 12 percent higher maximum wave heights in the near-field and slightly lower tsunamis in the far-field compared to a static sourcing, while the waves arrive later. The same scenario is further extended to a large-scale, fully coupled 3-D earthquake-tsunami model that captures dynamic rupture, seismic and acoustic wave propagation, and time-dependent tsunami generation. The fully coupled simulation reveals complex interactions between acoustic and tsunami waves during the early generation phase, including dispersion and wave conversions between seismic and ocean-acoustic waves, which static or linked models do not capture. Overall, our results provide quantitative evidence that source complexity, both spatial and temporal, affects tsunami near–field amplitudes, far–field decay, and arrival metrics, with important implications for physics–based tsunami hazard assessment along the Libyan and Egyptian coasts, the southern Aegean, and eastern Mediterranean shorelines. These results also highlight the value of integrating 3-D dynamic rupture modelling with tsunami simulations and enhance our understanding of tsunami generation mechanisms. 

Abstract Tsunami wave observations far from the coast remain challenging due to the logistics and cost of deploying and operating offshore instrumentation on a long-term basis with sufficient spatial coverage and density. Distributed Acoustic Sensing (DAS) on submarine fiber optic cables now enables real-time seafloor strain observations over distances exceeding 100 km at a relatively low cost. Here, we evaluate the potential contribution of DAS to tsunami warning by assessing theoretically the sensitivity required by a DAS instrument to record tsunami waves. Our analysis includes signals due to two effects induced by the hydrostatic pressure perturbations arising from tsunami waves: the Poisson’s effect of the submarine cable and the compliance effect of the seafloor. It also includes the effect of seafloor shear stresses and temperature transients induced by the horizontal fluid flow associated with tsunami waves. The analysis is supported by fully coupled 3-D physics-based simulations of earthquake dynamic rupture, seismo-acoustic waves, and tsunami (i.e., gravity wave) propagation. The strains from seismo-acoustic waves and static deformation near the earthquake source are orders of magnitude larger than the tsunami strain signal. We illustrate a data processing procedure to discern the tsunami signal. With enhanced low-frequency sensitivity on DAS interrogators (instrumental self-noise floor$$\approx 1 \times 10^{-12}$$ $\approx 1\times {10}^{-12}$ ), we find that, on seafloor cables located above or near the earthquake source area, tsunamis are expected to be directly observable with a sufficient signal-to-noise ratio within a few minutes of the earthquake onset. These results pave the way towards faster tsunami warning enabled by seafloor DAS. 

SUMMARY We simulate an instantaneous time mirror (ITM), that is, a rapid short-duration change in elastic material properties, using numerical experiments in time-varying isotropic elastic media. Our implementation in the seismic wave propagation software SeisSol is based on high-order discontinuous Galerkin discretization with ADER time stepping. We develop an eigenvector-based analytical solution for time interfaces for general linear hyperbolic wave systems and apply it to analyse the energy balance at time boundaries and ITMs. The energy increases for all intermittent medium changes for all impedance scaling factors. Our numerical implementation is validated against these analytical solutions and achieves high-order convergence. Its accuracy is further corroborated by estimates of reflection and transmission coefficients and observed frequency shifts across time boundaries, and by acoustic wave speed estimates obtained from focal spots associated with ITM-generated converging P waves that are consistent with theoretical predictions and ground truth values, respectively. We use the ITM implementation to simulate the partitioning of seismic body waves excited by a point source in a spatially homogeneous elastic full space. The response to an intermittent short change in the elastic parameters yields a diverging and converging P and S wavefield. A systematic scaling of the elastic parameters is then used to steer independent ITM reflections of either P or S waves. Numerical ITM solutions as developed here can be used to synthesize converging wavefields in seismic imaging applications, and more generally to analyse the behaviour and manipulation of seismic wavefields in space-time varying media.