Kevin R. Milner

We present an update to the 2023 U.S. 50‐state National Seismic Hazard Model, based on questions and comments received from a range of users, as well as additional internal and building code reviews since its initial release. For the conterminous United States, the update includes (1) revised weighting of ground motion models in the Intermountain West, (2) revised weighting of seismicity models across the stable crust of the central and eastern United States, (3) corrections to fault rupture input files and hazard codes, and (4) adjustments to stable crust nonlinear site amplification terms at soft soil sites. The model also updates the implementation of moment magnitude 7–8 subduction interface earthquake ruptures for Alaska and includes an alternate, full‐catalog rate model for Hawaii. The revised model is available through the U.S. Geological Survey online tools and is accompanied by a data release of high‐resolution hazard datasets for eight site conditions (from soft soil to hard rock with a time‐averaged shear‐wave velocity, V _S30 , ranging from 150 to 1500 m/s) in terms of peak horizontal ground acceleration, peak ground velocity, and 5% damped pseudo‐spectral accelerations at 21 spectral periods. It has been approved as the basis for updated design provisions and is being considered for inclusion in other hazard and risk assessment products.

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.

The Statewide California Earthquake Center has developed the CyberShake platform to perform 3D simulation-based probabilistic seismic hazard analysis. Using earthquake event rates, kinematic ruptures, and a 3D seismic velocity model, all constrained by observations, CyberShake utilizes wave propagation simulations and a reciprocity-based approach to generate hundreds of thousands of seismograms per site of interest. These simulation results are used to produce hazard data products including uniform hazard spectra, significant durations, site-specific hazard curves, and regional hazard maps. Previous CyberShake modeling efforts have computed seismic hazard models in southern California and Central California. In this article, we present CyberShake Study 24.8, which implements the CyberShake workflow for 315 sites in the greater San Francisco Bay area, extending physics-based hazard estimates to the entire San Andreas fault system. We construct the velocity model by tiling together three 3D and 1D models and applying a near-surface VS30-based merged taper to the top 700 m. We used the AWP—ODC wave propagation code to perform the ground-motion simulations, and the Graves–Pitarka kinematic rupture generator to create slip time histories for individual events. CyberShake Study 24.8 produces overall slightly reduced hazard compared to the NGA-West2 ground-motion models (GMMs). Specifically, CyberShake produces lower hazard immediately surrounding the San Francisco Bay likely due to use of a minimum VS of 400 m/s, lower hazard in near-fault regions, and higher hazard at long periods near sedimentary basins in the velocity model. We find larger directivity effects for CyberShake in northern California than in southern California, likely due to increased uniformity in fault characteristics and orientation in northern California. Validation results using the 1989 M 6.9 Loma Prieta event show that CyberShake reproduces observed ground motions with reduced bias compared with GMMs and 1D simulations.

The characteristic magnitude distribution hypothesis predicts a higher rate of large earthquakes than a Gutenberg–Richter extrapolation of the small-earthquake rate would imply. Characteristic magnitude distributions have been commonly applied to faults in probabilistic seismic hazard analysis (PSHA), and in modern models they can emerge from the way short-term seismicity constraints are combined with long-term geologic and geodetic constraints. We test the characteristic magnitude distribution hypothesis by comparing the fault-based magnitude distributions from the 2023 update to the National Seismic Hazard Model (NSHM23) in the Western United States with observed seismicity over the past 93 yr. We find that observed magnitude distributions fall outside the model-predicted confidence bounds in regions where NSHM23 produces characteristic magnitude distributions: in these regions, the model predicts higher rates of large earthquakes than are observed. An analysis of the earlier California model (Uniform California Earthquake Rupture Forecast, version 3) also reveals discrepancies between the modeled and observed magnitude distributions. In addition, we find that observed magnitude distributions near modeled faults are not significantly different from those in background regions. These results challenge the prevalence of characteristic magnitude distributions in fault-based seismic hazard models and call for a reassessment of how disparate data sets are integrated in PSHA.

A key use-case of geologic slip rates is within deformation models used in probabilistic seismic hazard analyses. Field-derived geologic slip rates have formed the cornerstone of deformation models in such applications for decades. Recent advancements in seismic hazard analyses have expanded the use of faults for which geologic slip rates are not well constrained using categorical slip rate estimates. Because of these advancements, application of a geologic deformation model for use in 2025 U.S. Geological Survey National Seismic Hazard Model Puerto Rico-U.S. Virgin Islands (NSHM25-PRVI) proved challenging due to: (1) a lack of field-based geologic slip rates, and (2) a lack of epistemic uncertainty distributions within a broad range of estimated slip rates. Preliminary versions of the NSHM25-PRVI model sampled these slip rate bins in a coincident manner along preferred and extreme value branches, which yielded untenable correlations in mean hazard results. To minimize the influence of correlated uncertainties amid these challenges, we develop a synthetic epistemic uncertainty distribution for deformation rate on each crustal fault. Each fault has a weighting schema across four possible distribution shapes: uniform, normal, triangular favoring local minima, and triangular favoring local maxima. The synthetic distributions are then sampled several times for each logic tree branch. The results provide a more realistic distribution of rates across the study region as compared with using correlated extrema sampling. This exploration of our method in a small region like PRVI can pave the way for larger-scale, more complicated applications (e.g., western United States).

We present the 2025 U.S. Geological Survey Puerto Rico and U.S. Virgin Islands (PRVI) time-independent earthquake rupture forecast (ERF), developed for the 2025 update to the National Seismic Hazard Model (NSHM) for PRVI. The updated ERF improves upon a prior model from 2003, including an expanded fault inventory with slip-rate estimates, updated seismicity catalogs, and refined subduction zone geometries and deformation models. It applies the fault-system inversion methodology to solve for rates of ruptures on modeled faults, adapted from the 2023 NSHM (NSHM23) for the western United States, including the first application of the inversion to model rates on a U.S. subduction interface. Off-fault and intraslab seismicity are constrained by observed seismicity and use updated methods developed for NSHM23. Uncertainties in model components are substantial, and the ERF represents epistemic uncertainties through a comprehensive logic tree consisting of 1.7 billion logic-tree branches combined across all sources.

Earthquake rupture forecasting is a critical component of seismic hazard analysis and requires identifying physically plausible rupture paths across complex fault networks. At its core, this task could be formulated as a large-scale combinatorial optimization problem, which involves selecting an optimal subset of fault segments from a set of candidates. Such problems pose significant challenges for traditional algorithms, as the number of admissible multi-fault ruptures grows combinatorially. Current operational workflows rely on locally applied plausibility filters. While computationally efficient, such greedy local heuristics risk excluding globally competitive rupture scenarios, particularly in regions with dense fault connectivity and competing rupture pathways. This work investigates whether quantum annealing hardware can serve as a sampling accelerator for exploring the ensemble of physically plausible ruptures beyond what is currently accessible to classical approaches. Designing a problem formulation that remains scientifically meaningful while respecting the constraints of current quantum hardware (size, noise, etc.) is nontrivial, and naïve encodings often collapse under these limitations. We encode rupture plausibility modeling as low energy solutions to a quadratic unconstrained binary optimization (QUBO) problem defined on a fault-network middle graph. Binary variables represent activated fault segments, while local interaction terms encode Coulomb stress transfer (as a proxy for pairwise rupture compatibility), continuity preferences, and branching behavior. Here, we adopt a classical-quantum hybrid workflow: the quantum annealer is used to sample globally competitive rupture candidates, while classical post-processing implements physical constraints to filter out non-physical ruptures. The workflow is demonstrated on a subset of the fault network used in the 2023 USGS National Seismic Hazard Model including over 200 fault segments, using both simulated thermal and quantum annealing on D-Wave hardware. Results show that our hardware-aware formulation and conditioning enable robust sampling on comparatively large fault-network instances. Views on quantum computing are polarized: some overstate its power, while others dismiss it as impractical. Our results help bridge this by establishing a concrete, testable pathway for integrating quantum annealing into rupture-modeling workflows on existing purpose-built quantum devices.—---This work is supported by the ICSC National Research Centre for High Performance Computing, Big Data and Quantum Computing (CN00000013, CUP D53C22001300005) within the European Union-NextGenerationEU program (National Recovery and Resilience Plan (PNRR) - Mission 4 Component 2 Investment 1.4.)

We articulate a scientific vision and roadmap for the development of improved Earthquake Rupture Forecast (ERF) models, which are one of the two main modeling components used in modern seismic hazard and risk analysis. One primary future objective is to provide fully time-dependent models that include both elastic rebound and spatiotemporal clustering nationwide, which can be particularly important for shorter-term hazard and risk considerations (e.g., earthquake insurance products). We also discuss the importance and perennial challenges associated with quantifying epistemic uncertainties, including those associated with deformation-model slip rates, unquantified sampling errors with respect to off-fault seismicity, and any spatial covariances. The need for more physics-based approaches is also emphasized, as is the benefit of adding model valuation (quantifying usefulness) to our verification and validation protocols. Given the multidisciplinary and system-level nature of this activity, modular design is critical. Future updates will also draw from best-available science by both the U.S. Geological Survey and the external community. The primary goal of this article is to highlight plans that guide research and facilitate community engagement with model development, especially with respect to lowering the entry barrier for early-career scientists and engineers. The article is written so readers can focus on the sections that interest them most (see Manuscript Organization section), with the Introduction and the Discussion section providing a stand-alone overview and summary, respectively.

We use the deterministic earthquake simulator RSQSim to generate complex sequences of ruptures on fault systems used for hazard assessment. We show that the source motions combined with a wave propagation code create surface ground motions that fall within the range of epistemic uncertainties for the Next Generation Attenuation-West2 set of empirical models. We show the model is well calibrated where there are good data constraints, and has good correspondence in regions with fewer data constraints. We show magnitude, distance, and mechanism dependence all arising naturally from the same underlying friction. The deterministic physics-based approach provides an opportunity for better understanding the physical origins of ground motions. For example, we find that reduced stress drops in shallow layers relative to constant stress drop with depth lead to peak ground velocities in the near field that better match empirical models. The simulators may also provide better extrapolations into regimes that are poorly empirically constrained by data because physics, rather than surface shaking data parameterizations, is underlying the extrapolations. Having shown the model is credible, we apply it to a problem where observations are lacking. We examine the case of crustal faults above a shallow subduction interface seen to break coseismically in simulations of the New Zealand fault system. These types of events were left out of consideration in the most recent New Zealand national seismic hazard model due to the modeling complexity and lack of observational data to constrain ground-motion models (GMMs). Here, we show that in the model, by breaking up the coseismic crustal and interface rupturing fault motions into two separate subevents, and then recombining the resulting ground-motion measures in a square-root-of-sum-of-squares incoherent manner, we reproduce well the ground-motion measures from the full event rupture. This provides a new method for extrapolating GMMs to more complex multifault ruptures.