Estimating the likelihoods of many earthquake scenario ruptures in a region in hazard models: Making grand inversions Bayesian

A robust seismic hazard model for a region, in principle, requires a sense of the likelihood of every conceivable earthquake affecting that region - or, short of that, the likelihood of an exhaustive (with respect to hazard) set of possible earthquake scenarios. A central component of this (as implemented in recent seismic hazard models for California, New Zealand, the United States and elsewhere) is a "grand inversion" approach (Page et al., 2014; Field et al., 2014, 2021; 2024; Milner et al., 2022; Milner and Field, 2024), in which one:

1) generates a large (order 1e5-1e6) set of simple scenario ruptures on known faults, noting their magnitudes and how much surface slip each would produce (model matrix G);

2) assembles geologic and geodetic constraints on total fault-slip rates, paleoseismologic constraints on large-earthquake recurrence intervals, and seismic-catalogue constraints on the total magnitude-frequency distribution in the system (constraint vector d);

3) "inverts" the constraint vector and model matrix to estimate the rate of each scenario rupture in the system.

The model space is large and underdetermined (Page et al., 2014), so up to now, a simulated-annealing approach has been used to efficiently find a global-minimum solution that best fits the constraints. Then the uncertainties on the constraints, trade-offs between model elements, and prediction uncertainties have been propagated into the solution space by carrying many grand inversions with different input constraints (e.g. geologic or geodetic data or both, various b-values, various slip scaling laws) in each branch of a large logic tree, and by toggling importance weights on different constraints.

These grand inversions formed a central component of the New Zealand National Seismic Hazard Model 2022 (Gerstenberger et al., 2024). In this process, we identified two characteristics that merit further work and may substantially impact estimated hazard levels. First, the current simulated-annealing approach return very sparse solutions. The input scenario-rupture sets for New Zealand feature several hundred faults, several thousand fault subsections and 1e5-1e6 plausible ruptures, but the grand inversions typically assign a rate of 0 to all but ~1000 ruptures, and many faults have only one or a few nonzero-rate ruptures (effectively a nearly characteristic model). Second, the grand inversions output only the global-minimum solution rather than the entire model-space exploration. This means that many of the uncertainties on the constraints (those that are not toggled overhead as alternate logic-tree branches), such as fault slip rate uncertainties, are not propagated into the model space except as weights. To overcome these limitations, we are making the grand inversions Bayesian by replacing the simulated-annealing approach with a Monte Carlo search with No U-Turn Sampling, and outputting the full posterior distribution of the model space. This will allow the grand inversions to propagate the full range of possible scenario ruptures on each fault into hazard estimates (modulo the data constraints) rather than only a few select ruptures.