# Instantaneous physics-based ground motion maps using reduced-order modeling

###### Abstract

Physics-based simulations of earthquake ground motion are useful to complement
recorded ground motions. However, the computational expense of performing
numerical simulations hinders their applicability to tasks that require real-
time solutions or ensembles of solutions for different earthquake sources. To
enable rapid physics-based solutions, we present a reduced-order modeling
(ROM) approach based on interpolated proper orthogonal decomposition (POD) to
predict peak ground velocities (PGVs). As a demonstrator, we consider PGVs
from regional 3D wave propagation simulations at the location of the 2008 Mw
5.4 Chino Hills earthquake using double-couple sources with varying depth and
focal mechanisms. These simulations resolve frequencies $\leq$ 1.0 Hz and
include topography, viscoelastic attenuation, and S-wave speeds $\geq$ 500
m/s. We evaluate the accuracy of the interpolated POD ROM as a function of the
approximation method. Comparing the radial basis function (RBF), multilayer
perceptron neural network, random forest, and $k$-nearest neighbor, we find
that the RBF interpolation gives the lowest error ($\approx$ 0.1 cm/s) when
tested against an independent dataset. We also find that evaluating the ROM is
$10^{7}-10^{8}$ times faster than the wave propagation simulations. We use the
ROM to generate PGV maps for one million different focal mechanisms, in which
we identify potentially damaging ground motions and quantify correlations
between focal mechanism, depth, and accuracy of the predicted PGV. Our results
demonstrate that the ROM can rapidly and accurately approximate the PGV from
wave propagation simulations with variable source properties, topography, and
complex subsurface structure.

Journal of Geophysical Research: Solid Earth

Scripps Institution of Oceanography, University of California San Diego, La
Jolla, USA Department of Earth and Environmental Sciences, Ludwig-
Maximillians-Universität München, Munich, Germany John M.
<EMAIL_ADDRESS>
Physics-based simulations are useful for seismic hazard assessment and ground
motion prediction but are computationally expensive.

Our reduced-order model accurately predicts simulated peak ground velocities
for scenarios, including 3D wavefield and topographic effects.

The reduced-order model produces physics-based ground motion maps in
milliseconds for variable source depths and focal mechanisms.

## Plain Language Summary

Computer simulations can be used to predict the intensity of ground shaking
caused by earthquakes. However, these simulations require significant amounts
of computing time to obtain accurate results, making them too slow for real-
time systems. In this work, we apply a technique that allows us to obtain
instantaneous and accurate approximations of a scenario earthquake simulation
using simplified models. The simplified models are fast enough to evaluate
this earthquake scenario on demand. We use the simplified models to predict
shaking intensities for one million different fault planes at the location of
the 2008 Chino Hills earthquake. By analyzing the predicted shaking
intensities from these earthquakes, we identify possible areas of strong
shaking, areas that might otherwise be unidentified if fewer earthquake
scenarios were examined. These results suggest that these types of simplified
models could be considered to assess regional earthquake hazards and to inform
shaking intensity predictions in earthquake early warning.

##  1 Introduction

A fundamental challenge in seismology and earthquake engineering is the
estimation of the ground shaking caused by earthquakes, especially for areas
with high population densities or containing critical infrastructure and in
the near-source region where observations are sparse. For example, after large
earthquakes, it is important to rapidly characterize the shaking intensity to
identify areas of potential damage and guide rapid response efforts.
Accurately estimating shaking intensities in near-real-time is one of the
primary goals of the ShakeMap product developed by the U.S. Geological Survey
[Worden Wald (2016)]. Additionally, real-time estimates of the expected ground
motion intensity are needed in earthquake early warning (EEW) to issue alerts
to people and critical systems, ideally in advance of the strongest shaking
[Kohler . (2017)]. Also, systematically quantifying the potential level of
ground shaking from as of yet unobserved earthquakes that may someday occur is
important for earthquake preparedness and seismic hazard assessment
[<]e.g.,¿[]gravesCyberShakePhysicsbasedSeismic2011,Baker2014,bommerSelectionGroundmotionPrediction2010.

Several strategies have been developed to estimate earthquake ground motions.
Instrumental recordings and non-instrumental data such as “Did You Feel It?”
reports [Wald . (2011)] directly inform the amplitudes and spatial
distribution of seismic ground motion. However, the generally sparse spatial
coverage of observations leads to significant uncertainties in the shaking
distributions [Wald . (2008)]. To fill in the gaps, recent approaches include
spatial and spectral interpolation [Worden . (2018)], Gaussian process
regression [Tamhidi . (2021)], and conditional simulation [Engler . (2022),
Bailey . (2022)], but rapidly estimating the full distribution of shaking
intensity remains a challenge. Most EEW systems and large-scale seismic hazard
studies
[<]e.g.,¿[]stirlingNationalSeismicHazard2012,kohlerEarthquakeEarlyWarning2017,petersen2018UpdateUS2020,
melettiNewItalianSeismic2021 instead use empirical ground motion models (GMMs)
which give instantaneous ground motion estimates for given source, path, and
site parameters. GMMs rely on the availability of a large ground motion
dataset recorded in a specific region or tectonic setting
[<]e.g.,¿[]bozorgniaNGAWest2ResearchProject2014,bozorgniaNGASubductionResearchProgram2021,gouletNGAEastGroundMotionCharacterization2021,
which may not always be available. Thus, GMMs typically rely on the ergodic
assumption [Anderson Brune (1999)] where all the data from a certain tectonic
environment are grouped together to develop a single model that describes the
average scaling with magnitude, distance, and site conditions. However, the
uncertainty of GMM predictions can be reduced with nonergodic source, path, or
site correction terms if they are well-constrained
[<]e.g.,¿[]baltayUncertaintyVariabilityEarthquake2017,sahakianGroundMotionResiduals2019,parkerEmpiricalMapBasedNonergodic2022.
Going beyond the source-based approach of GMMs, other methods have been
developed that directly use ground motion or wavefield observations to
forecast the shaking at distant locations [Hoshiba Aoki (2015), Kodera .
(2018), Furumura . (2019), Cochran . (2019), Oba . (2020)]. A recent
wavefield-based approach developed by nagata2023seismic uses reduced-order
modeling to reconstruct seismic wavefields from sparse measurements, although
their approach does not include the source parameterization that is needed to
solve inverse problems.

Complementary to GMM and wavefield-based prediction techniques, numerical
simulations of earthquake rupture and seismic wave propagation provide
physics-informed estimates of shaking intensity
[<]e.g.,¿[]aagaardCharacterizationNearsourceGround2001,
olsenStrongShakingAngeles2006,moschetti3DGroundmotionSimulations2017,
rodgersRegionalScale3D2020,taufiqurrahmanBroadbandDynamicRupture2022a. These
simulations can assume a kinematic description of the earthquake source
[<]e.g.,¿[]gravesKinematicGroundmotionSimulations2016, or solve for the
earthquake slip history with spontaneous dynamic rupture models
[<]e.g.,¿[]rippergerVariabilityNearfieldGround2008,withersValidationDeterministicBroadband2018.
In both cases, the seismic wave equation is solved to simulate the propagation
of seismic waves within the Earth. Physics-based 3D ground motion simulations
have advanced in terms of realism and agreement with real data due to
increasingly dense observations, a better understanding of Earth structure and
ever-increasing computing power. For example, recent advances in community
velocity model (CVM) development in southern California have been used to
inform physics-based seismic hazard assessment [Graves . (2011), Milner .
(2021)]. Important progress has been made in generating broadband synthetics
which usually involve a stochastic component for the high frequencies
[<]e.g.,¿[]booreStochasticSimulationHighfrequency1983,frankelBroadbandSyntheticSeismograms2018.

Deterministically resolving the high-frequency wavefield makes physics-based
ground motion modeling computationally expensive, typically requiring high-
performance computing (HPC) technology
[<]e.g.,¿[]Cui2013,Heinecke2014,rodgersRegionalScale3D2020,Savran2020,pitarkaDeterministic3DGroundMotion2022b.
The high cost of each high-resolution numerical simulation challenges
applications in which many forward models must be evaluated, such as dynamic
earthquake source inversion problems
[<]e.g.,¿[]gallovicBayesianDynamicFinitefault2019b or physics-based
probabilistic seismic hazard assessment. Even though modeling simplifications
help reduce the computational expense (e.g., by setting an upper threshold for
the minimum S-wave velocity, $V_{S}$), regional-scale broadband simulations
can require up to $10^{3}-10^{4}$ CPU hours (CPUh) of computing time
[<]e.g.,¿[]huHzDeterministic3D2022,taufiqurrahmanBroadbandDynamicRupture2022a.

To address these challenges, recent efforts have enabled the rapid generation
of synthetic seismograms using Green’s functions or deep learning. For global
seismology, 1.0 Hz discrete Green’s functions for axisymmetric Earth models
are available on-demand via InstaSeis [van Driel . (2015)], and 3D Green’s
functions based on the spectral element method [Komatitsch Tromp (1999)] are
typically available a few hours after $M_{w}>5.5$ earthquakes [Tromp .
(2010)]. However, the limited resolution of global approaches are not suitable
for regional ground motion applications. At the regional scale,
wangMoving1D3D2020 created a 0.5 Hz 3D Green’s function database for stations
in the southern California seismic network to enable real-time moment tensor
inversions with 3D models. However, such Green’s function databases often
require several terabytes of storage, presenting a major challenge in
accessibility. Other approaches used deep learning to learn approximate
solutions to the wave equation and produce seismograms for arbitrary velocity
models [Moseley . (2018), Moseley . (2020), Smith . (2020), Yang . (2021),
Yang . (2023), Rasht-Behesht . (2022)] or create neural network GMMs based on
simulations [Withers . (2020)]. However, machine-learning approaches have not
yet directly reproduced 3D path or topographic effects that are important for
ground motion prediction.

Here, we present a surrogate modeling approach that aims to address the need
for rapid generation of 3D physics-based synthetics and ground motion
intensity maps. Our approach is based on a reduced-order model (ROM), a
particular type of surrogate model. In the model order reduction framework, a
high-fidelity forward model can be approximated by a surrogate model, which
typically has many fewer degrees of freedom and is fast to evaluate. As a
result, a well-trained ROM can enable real-time and many-query evaluation of a
system and can closely approximate the solution of an expensive forward
simulation. We predict the peak-ground velocity (PGV) from a scenario of
regional ground motion simulation. Our method produces physics-informed
estimates instantaneously, thus, at the same speed as GMMs. While ROMs have
rarely been applied in geophysics thus far, they have been demonstrated to
accurately approximate solutions to linear and non-linear physics-based models
in other fields of science and engineering
[<]e.g.,¿[]willcoxBalancedModelReduction2002,bennerSurveyProjectionbasedModel2015,hesthavenReducedBasisMethods2022.

The structure of this paper is organized as follows. We first provide an
overview of our ROM development workflow based on an interpolated proper
orthogonal decomposition (iPOD) technique in Section 2.1. We use simulations
to create simulated PGV maps for different earthquake sources and decompose
them using iPOD. We create the ROM by fitting function approximators that
predict the iPOD coefficients for different values of the input source
parameters. We can then use the ROM to instantaneously predict PGV maps for
new earthquakes that were not modeled. We describe the details of the forward
problems for both a 1D layered and a 3D velocity model with topography in
Section 2.2. We describe how we will use four different function approximators
to predict the simulation output in Section 2.3. In our results Section 3, we
evaluate the accuracy of the ROMs using an independent testing dataset. We
provide illustrative examples highlighting the ROM’s ability to reproduce PGV
maps faithfully. Finally, we present a potential application and sensitivity
analysis emphasizing the importance of hypocentral depth in PGV predictions.
Lastly, we discuss in Section 4 how reduced-order modeling may path new
avenues for physics-based hazard analysis and for using high-performance
computing empowered simulations in real-time systems by enabling the
instantaneous generation of physics-informed ground motion maps.

##  2 Methods

###  2.1 Reduced-order modeling workflow

We consider the earthquake source to be a point source in space, parameterized
by its hypocentral depth, strike, dip, and rake. These four parameters will be
denoted by $\bm{p}$. The ground motion intensity metric we use is the surface
peak-ground velocity (PGV) denoted by $\bm{q}$. We define the mapping between
the parameters $\bm{p}$ and the PGV $\bm{q}$ using a physics-based approach
given by the (visco-)elastic wave equation for an isotropic medium. When a 3D
velocity model is used to describe the Earth’s crustal structure and/or non-
planar topography is considered, this problem must be solved numerically. The
computational cost of performing these numerical wave propagation simulations
precludes their direct usage within real-time or rapid response
infrastructure, as well as many-query tasks, for example as required for
global sensitivity analysis and forward uncertainty quantification.

Referring to the numerical method used to solve the wave equation and produce
the PGV $\bm{q}$ as the full-order model (FOM), we wish to construct a
reduced-order model (ROM) that approximates $\bm{q}$ given the earthquake
parameters $\bm{p}$ as input (Figure 1). While solving the FOM is a
computationally expensive problem that requires HPC, the ROM has fewer degrees
of freedom resulting in lower computational complexity compared to the FOM.
Thus, the ROM is much faster to evaluate than the FOM. To inform the ROM, we
obtain FOM solutions for a set of points $\\{\bm{p}_{1}$, $\bm{p}_{2}$,
$\ldots$, $\bm{p}_{n_{s}}\\}$ where $n_{s}$ is the total number of simulations
that will be performed. Each parameter vector $\bm{p}_{i}$ is contained within
the four-dimensional space $\mathcal{P}$ which encompasses the range of source
depths and focal mechanism angles we consider; the ROM is only applicable for
sources contained within $\mathcal{P}$. We then split the set of FOM solutions
into two groups: training FOM solutions and testing FOM solutions. After using
the training FOM solutions to train the ROM, we validate it against the
testing FOM solutions. If there is a good agreement, then the ROM can be used
as an efficient substitute, or surrogate, for the FOM.

Figure 1: Conceptual illustration of the reduced-order model (ROM) development
and validation procedures. We first obtain full-order model (FOM) solutions
for points in the parameter space $\mathcal{P}$ and divide them into training
(red squares) and testing solutions (blue squares). Here, FOM refers to the
forward problem of a numerical solution of the wave equation and producing a
surface peak ground velocity (PGV) map. After training the ROM using the
training FOM solutions, one can evaluate the ROM for any new point in
$\mathcal{P}$ (open circles). The ROM solutions are validated against the
testing FOM solutions (overlapping circles and blue squares). As shown on the
right, the FOM and ROM output similar PGV maps that depend on the input source
parameters; the ROM, however, is a computationally inexpensive model with
fewer degrees of freedom.

###  2.2 Forward model

Our forward problem is the simulation of PGV maps that result from earthquakes
with a given hypocentral depth and focal mechanism. We numerically simulate
the 3D propagation of seismic waves produced by each source, varying the
strike, dip, rake, and hypocentral depth of the source in each simulation. Our
source parameterization is motivated by seismicity in Southern California
[SCEDC (2013), Field . (2014)] but simplified for the proof-of-concept scope
of this study. As a source-time function that describes the rate of earthquake
moment release, we use $\dot{M}(t)=\frac{M_{0}t}{T^{2}}e^{-t/T}$ with $T=0.4$
s. This results in a Brune-type source spectrum with a stress drop of
approximately 5.0 MPa [Brune (1970)]. We do not vary the total scalar seismic
moment, prescribing $M_{w}=5.4$ in all simulations. We assume point-source,
double-couple moment tensor earthquake source descriptions, which is a
reasonable assumption for most small earthquakes resulting from shear faulting
[Julian . (1998)]. The hypocentral depths vary from 2.0 to 20.0 km.

To generate ground motion maps from the seismic wave simulations for each
source, we first define a high-resolution area of interest. We select a 30 km
by 30 km square area (Figure 2a) centered at the epicenter of the 2008 $M_{w}$
5.4 Chino Hills earthquake (33.949∘ N, 117.776∘ W), an area that has
potentially significant topographic and basin effects on ground motion
[Taborda Bielak (2013)]. Note that the complete computational mesh is larger
(150 $\times$ 150 $\times$ 100 km) than the area of interest to prevent
artificial boundary reflections from affecting the PGV measurements. We limit
the computational expense of the forward simulations by selecting
$f_{max}=1.0$ Hz as the maximum frequency to fully resolve.

For each earthquake source, we perform forward simulations with three
different setups of increasing complexity (Figure 2b). The first setup uses a
flat, layer-over-halfspace velocity model, which defines a 1D velocity
structure (denoted as “LOH”, see Table S1 for the material parameters). The
second setup “3D-1400” uses a 3D community velocity model for Southern
California, SCEC CVM-H v15.1.0 [Shaw . (2015)], and includes topography.
“1400” indicates that we set a floor for the minimum S-wave velocity at
$V_{S,\text{min}}=1400$ m/s. In both LOH and 3D-1400, we assume linear,
elastic wave propagation. The third setup “3D-500A” is similar to 3D-1400
except that $V_{S,\text{min}}=500$ m/s and that we account for frequency-
independent viscoelastic attenuation with $Q_{p}=0.1V_{S}$ and
$Q_{s}=0.05V_{S}$ (for $V_{S}$ measured in m/s) [K. Olsen . (2003)]. The
viscoelastic attenuation setup assumes a Maxwell material with three damping
mechanisms with their respective relaxation frequencies of 0.05, 0.5, and 5.0
Hz [Uphoff . (20171)]. These different modeling setups allow us to examine how
the ROM predictions and accuracies are affected by the modeling assumptions
related to each simulation regime.

Figure 2: (a) Map of the area of interest (red box) and the surrounding region
in Southern California. The epicenter of the simulated earthquakes (which is
the same as the 2008 $M_{w}$ 5.4 Chino Hills earthquake) is indicated by a
black star and a site of interest is indicated by a black triangle. (b) Cross-
sections (A - A’) of the velocity models used for wave propagation simulation
of the three forward models (LOH, 3D-1400, and 3D-500A). Note that the cross-
sections correspond to the area of interest, though our simulations use larger
meshes to prevent any potential spurious boundary reflections. As shown in the
cross-sections, we vary the source depth and focal mechanism for each of the
$n_{s}=5000$ simulations. The source depths range between $z_{min}$=2.0 km and
$z_{max}$=20.0 km, though the cross-sections show a maximum depth of 5.0 km to
better show the near-surface velocity structures.

To perform the seismic wave propagation simulations for each source, we use
SeisSol (https://github.com/SeisSol), an open-source code based on the
arbitrary high-order accurate derivative discontinuous Galerkin (ADER-DG)
method [Dumbser Käser (2006)] with end-to-end optimization for HPC
infrastructure [<]e.g.,¿[]Rettenberger2016,Uphoff2017,Krenz2021. SeisSol uses
fully non-uniform unstructured tetrahedral meshes that can statically adapt to
complex 3D subsurface structures and topography. We create 3D tetrahedral
meshes for the simulations using velocity-aware meshing, i.e., automatically
adapting SeisSol’s mesh resolution to the 1D or 3D velocity model (Figure 2b).
We ensure sufficient resolution using four elements with basis functions of
polynomial degree of 4 per shortest wavelength [Käser . (2008)]. To inform the
3D models, we obtain the CVM-H velocity model using the SCEC UCVM software
[Small . (2017)]. We obtain the topography data from the 1-meter resolution
digital elevation model from the U.S. Geological Survey 3D Elevation Program
(see data availability statement). We create the meshes with Gmsh [Geuzaine
Remacle (2009)], resulting in 550,001 tetrahedral elements for LOH, 129,888
elements for 3D-1400, and 1,189,840 elements for 3D-500A. All simulations were
performed using a degree four polynomial in space and time (translating to a
fifth-order accurate scheme in space and time) and were executed on SuperMUC-
NG using 17 compute nodes. Each compute node consists of 48 cores employing
Intel Skylake Xeon Platinum 8174 processors. Using these resources the
runtimes per simulation were 59.5 s, 23.0 s, and 390.0 s for the LOH, 3D-1400,
and 3D-500A models, respectively. Assuming perfect utilization of the
17$\times 48$ cores, these runtimes translate to 13.5, 5.2, and 88.4 CPUh. We
save the three-component free-surface velocity time-series for $n_{r}=3600$
receivers with a spacing of 500 m and sampling rate of 10.0 Hz. To create a
ground motion map for each earthquake source, the last step is to compute the
PGV from the time series of the horizontal components. We demean, detrend, and
apply a 1.0 Hz lowpass filter to the time series before combining the two
horizontal components using the orientation-independent RotD50 measure [Boore
(2010)].

###  2.3 Reduced-order model

We now describe in detail how to define the ROM that approximates the mapping
between the parameters $\bm{p}$ and the PGV map $\bm{q}$.

####  2.3.1 Parameter space sampling

To inform the ROM, we first select the points in the parameter space
$\mathcal{P}$ where we will obtain the FOM solutions. To determine these
points, we use a low-discrepancy, four-dimensional Halton sequence [Halton
(1960)]. The low-discrepancy property of this sequence means that it closely
approximates points drawn from a uniform distribution, giving us points that
evenly cover our entire parameter space. We use the first $n_{s}=5000$ points
from this sequence, where $n_{s}$ is the total number of earthquake sources
for which we evaluate the FOM. We gather the parameters determined by the
Halton sequence into a matrix
$\bm{P}=[\textbf{{p}}_{1},\textbf{{p}}_{2},\ldots,\textbf{{p}}_{n_{s}}]$ and
the PGV data generated by the FOM into a matrix
$\bm{Q}=[\textbf{{q}}_{1},\textbf{{q}}_{2},\ldots,\textbf{{q}}_{n_{s}}]$.
Generating the simulation data matrix $\bm{Q}$ for $n_{s}=5000$ sources
requires a total of 67.5 kCPUh, 26 kCPUh, and 442 kCPUh for the LOH, 3D-1400
and 3D-500A forward models, respectively. To create independent datasets for
training and testing of the ROM, we randomly separate the PGV maps (stacked in
$\bm{Q}$) that result from these sources into training
$(\bm{Q}_{\text{train}},\bm{P}_{\text{train}})$ and testing
$(\bm{Q}_{\text{test}},\bm{P}_{\text{test}})$ datasets, where the testing
dataset contains 10% of the simulated PGV maps ($n_{\text{train}}=4500$,
$n_{\text{test}}=500)$. The accuracy of the trained ROM can then be evaluated
using PGV maps from the independent testing dataset.

As the goal of the ROM is to interpolate the known PGV maps for new sources in
$\mathcal{P}$ that we did not simulate, we measure the Euclidean distances
between each testing parameter $\bm{p}\in\bm{P}_{\text{test}}$ and its
nearest-neighboring training parameter
$\bm{p_{\text{nearest}}}\in\bm{P}_{\text{train}}$, i.e.,

$d_{\text{nearest}}(\bm{p})=||\bm{p}-\bm{p_{\text{nearest}}}||_{2}.$

These are the distances after normalizing each dimension of the parameter
space (depth, strike, dip, and rake) to range between 0 and 1. In our testing
dataset, the distribution of distances has a mean value of 0.09 (Figure S1).
We also create a set of sources $\bm{P}_{\text{uniform}}$ that contains
$10^{6}$ uniformly distributed sources that evenly cover $\mathcal{P}$. While
it is not possible to evaluate the FOM to obtain error estimates for the
$10^{6}$ points in $\bm{P}_{\text{uniform}}$, we show that our testing
distribution is a good proxy of the errors as indicated by the mean and
standard deviations of the distances for both $\bm{P}_{\text{test}}$ and
$\bm{P}_{\text{uniform}}$ (Figure S1). For any earthquake source contained in
$\mathcal{P}$, the distance to the nearest training point will be, on average,
0.08. The maximum possible distance to any training FOM location is 0.17 units
which is our reported “fill distance” [De Marchi Schaback (2010)] for our
sampled FOM locations.

####  2.3.2 Interpolated proper orthogonal decomposition (iPOD)

We use a method for defining a ROM called the interpolated proper orthogonal
decomposition (iPOD) [<]e.g.,¿[]bui2003proper, druaultUseProperOrthogonal2005.
In this technique (Figure 3a), we first compute a modal decomposition for the
collected data using the singular value decomposition (SVD),

$\bm{Q}=\bm{U\Sigma V^{T}},$ (1)

where $\textbf{{U}}=[\bm{u}_{1},\dots,\bm{u}_{n_{s}}]$ contains the left
singular vectors (or modes) $\bm{u}_{i}$,
$\bm{\Sigma}=\text{diag}[\sigma_{1},\dots\sigma_{n_{s}}]$ is a diagonal matrix
containing the singular values $\sigma_{i}$, and V contains the right singular
vectors [Berkooz . (1993)]. We then define the matrix of POD coefficients
$\bm{A}$ to be,

$\displaystyle\bm{A}^{T}$ $\displaystyle=\bm{\Sigma V^{T}}$ (2)
$\displaystyle\bm{A}$ $\displaystyle=\bm{V}\bm{\Sigma},$ (3)

where we denote individual entries of $\bm{A}$ as $\alpha_{ij}$. Using Eqs.
(1) and (3) we can express every $\bm{q}_{i}\in\bm{Q}$ using a weighted sum
(Figure 3a) of the left singular vectors, i.e.,

$\bm{q}(\bm{p}_{i})=\bm{q}_{i}=\sum_{j=1}^{n_{s}}\alpha_{ij}\bm{u}_{j}.$ (4)

Eq. (4) is a restatement of the factorization defined in Eq. (1) and thus can
only be used to reconstruct values of $\bm{q}$ for parameters
$\bm{p}\in\bm{P}$. To enable predictions for parameters $\bm{p}\notin\bm{P}$,
the iPOD method utilizes a generalization of Eq. (4)

$\bm{q}(\bm{p})\approx\tilde{\bm{q}}(\bm{p})=\sum_{j=1}^{n_{s}}f_{j}(\bm{p})\bm{u}_{j},$
(5)

where each function $f_{j}(\bm{p})$, $j=1,\dots,n_{s}$ predicts the POD
coefficient for any parameter $\bm{p}$ associated with the $j^{\text{th}}$
mode (Figure 3b). The types of functions $f_{j}(\cdot)$ that we consider are
described in Section 2.3. After fitting the functions $f_{j}(\bm{p})$ with the
POD coefficients given by $\alpha_{ij}$, the PGV map $\tilde{\bm{q}}$ is
obtained by evaluating Eq. (5). In essence, iPOD builds an orthogonal basis
$\bm{u}_{i}$ based on simulation data $(\bm{Q})$ with which it will use to
define a predicted PGV map. The basis $\bm{u}_{i}$ are combined linearly, with
the associated weight of each basis defined by a functional approximation of
the coefficients $\bm{A}$. The POD coefficients needed to reconstruct the
training dataset are then given by $\bm{A}_{\text{train}}$.

Figure 3: Conceptual illustration of a reduced-order model based on
interpolated proper orthogonal decomposition (iPOD). We first decompose the
data matrix of PGV maps $\bm{Q}$ using the singular value decomposition (a).
The goal is then to fit a series of functions $f_{j}$ (b) that predict the
proper orthogonal decomposition (POD) coefficients $\alpha_{ij}$ for different
values of the input source parameters (e.g., hypocentral depth). The ROM-
predicted PGV maps are obtained by evaluating the functions to obtain
predicted POD coefficients, and then performing a weighted sum of the left
singular vectors $\bm{u}_{j}$.

####  2.3.3 Function approximation for $\alpha_{ij}$

We consider four types of function approximators to model the relationships
between the source parameters $\bm{p}$ and POD coefficients $\alpha_{ij}$.
These include radial basis functions, neural networks, random forests, and
$k$-nearest neighbors. We summarize these function approximation techniques
below.

  1. 1.

Radial basis function (RBF) interpolation is commonly used as an efficient
technique for interpolating high-dimensional, unstructured data [Lazzaro
Montefusco (2002)] and has been shown to be successful in other studies that
use iPOD
[<]e.g.,¿[]audouzeReducedorderModelingParameterized2009,xiaoNonintrusiveReducedorderModelling2015.
For RBF interpolation, the interpolator for the $j$-th mode $f_{j}(\cdot)$ is
given by

$f_{j}(\bm{p})=\sum_{i=1}^{n_{\text{train}}}w_{ij}\varphi(||\bm{p}-\bm{p}_{i}||)+\sum_{k=1}^{|V|}b_{kj}\psi_{k}(\bm{p})$
(6)

where the function $\varphi(\cdot)$ is an RBF, $w_{ij}$ in an entry of the
kernel weight matrix $\bm{W}$, $b_{kj}$ is an entry of the polynomial
coefficient matrix $\bm{B}$, and $\psi_{k}(\cdot)$ are monomial terms up to a
specified degree $d$ contained in the set $V$ with

$V=\Biggl{\\{}p_{1}^{q}p_{2}^{r}p_{3}^{s}p_{4}^{t}:\begin{cases}q+r+s+t\leq
d\\\ 0\leq q,r,s,t\leq d\end{cases}\Biggr{\\}}$ (7)

where $p_{1}$ denotes the hypocentral depth, $p_{2}$ denotes the strike, etc.
We determine the coefficients contained in $\bm{W}$ and $\bm{B}$ by solving
the linear system of equations

$\displaystyle\bm{\Phi W}+\bm{\Psi B}$ $\displaystyle=\bm{A}_{\text{train}}$
$\displaystyle\bm{\Psi}^{T}\bm{B}$ $\displaystyle=\bm{0}$

where

$\bm{\Phi}=\begin{bmatrix}\varphi(\left\lVert{\bm{p}_{1}-\bm{p}_{1}}\right\rVert)&\varphi(\left\lVert{\bm{p}_{2}-\bm{p}_{1}}\right\rVert)&\ldots&\varphi(\left\lVert{\bm{p}_{n_{\text{train}}}-\bm{p}_{1}}\right\rVert)\\\
\varphi(\left\lVert{\bm{p}_{1}-\bm{p}_{2}}\right\rVert)&\varphi(\left\lVert{\bm{p}_{2}-\bm{p}_{2}}\right\rVert)&\ldots&\varphi(\left\lVert{\bm{p}_{n_{\text{train}}}-\bm{p}_{2}}\right\rVert)\\\
\vdots&\vdots&\ddots&\vdots\\\
\varphi(\left\lVert{\bm{p}_{1}-\bm{p}_{n_{\text{train}}}}\right\rVert)&\varphi(\left\lVert{\bm{p}_{2}-\bm{p}_{n_{\text{train}}}}\right\rVert)&\ldots&\varphi(\left\lVert{\bm{p}_{n_{\text{train}}}-\bm{p}_{n_{\text{train}}}}\right\rVert)\\\
\end{bmatrix}$ (8)

and $\bm{\Psi}$ is the matrix of polynomials up to degree $d$. For the RBFs
$\varphi(\cdot)$, we test three types of kernels: thin plate spline
($\varphi(r)=r^{2}\ln(r)$), cubic ($\varphi(r)=r^{3}$) and quintic
($\varphi(r)=-r^{5}$). For each kernel, we solve the linear system using the
minimum polynomial degree (equal to one for the thin plate spline and cubic
kernels, and degree two for the quintic kernel). Note that we do not apply
smoothing, which implies that the POD coefficients are exactly reproduced for
the training dataset, i.e., for any $\bm{p}\in\bm{P}_{\text{train}}$. The type
of RBF kernel (i.e., thin plate spline, cubic or quintic) is a hyperparameter
that we tune using cross-validation.

  2. 2.

Artificial neural networks are another type of function approximator for the
POD coefficients [Hinton (1990), Hesthaven Ubbiali (2018)]. We use a fully
connected, feed-forward structure in which the neural network consists of an
input layer, two hidden layers, and an output layer. This type of structure is
also known as a multilayer perceptron (MLP). In this type of model, the input
layer receives the four earthquake source parameters from $\bm{p}$, and the
output layer returns the predicted POD coefficients. Each neuron in the hidden
layers receives an input from each neuron in the previous layer, evaluates an
activation function (ReLU), then passes the information forward through the
network to every neuron in the next layer using determined weights and biases.
We optimize the weights and biases by training the neural networks with the
Adam optimizer [Kingma Ba (2017)] and the squared error loss function.
Furthermore, we follow the method of swischukProjectionbasedModelReduction2019
to predict a subset of the POD coefficients which correspond to the $r$
largest singular values, choosing to train the neural networks using $r=$ 10,
30, and 50 modes from the SVD (Eq. (1)). The hyperparameters for MLP consist
of the rank $r$ and the number of neurons in each hidden layer.

  3. 3.

Random forests (RF) are a type of ensemble estimator in which many decision
trees are constructed, and the final estimate is the mean of each decision
tree’s individual estimate [Breiman (2001)]. We here use both 50 and 100
decision trees. The hyperparameters are given by the number of decision trees.

  4. 4.

$k$-nearest neighbors (kNN) is a relatively simple type of approximator in
which the predicted value is given as the mean of the $k$-nearest neighboring
points. We apply uniform weighting to each point in the neighborhood and use
the $\ell_{2}$ norm as the distance metric. We test using 3, 5, and 7
neighbors which define the hyperparameters.

####  2.3.4 Training procedure

Before training the function approximators RBF, MLP, RF, and kNN, we
standardize both the inputs (earthquake source parameters, $\bm{P}$) and
outputs (POD coefficients, $\bm{A}$) such that they have a mean of zero and
standard deviation of one. Note that we use all available modes (i.e., all
columns of $\bm{A}$) to train the ROMs, except for the MLP approximator where
we allow the number of modes to be a hyperparameter. To determine the values
of the hyperparameters, we use $k$-fold cross-validation with $k=5$. In this
procedure, the training dataset is further subdivided into 5 folds. For each
of the folds, the remaining $k-1$ folds are used to train the approximators,
and the fold that was not included is the validation dataset. To measure the
error of a ROM with a given set of hyperparameters, we use both the mean
absolute error and the mean absolute percentage error. We define the local (in
parameter space) mean absolute error as

${MAE}(\bm{p})=\frac{1}{n_{r}}\left\lVert{\bm{q}(\bm{p})-\tilde{\bm{q}}(\bm{p})}\right\rVert_{1}$
(9)

and the local mean absolute percentage error as

${MAPE}(\bm{p})=\frac{1}{n_{r}}\left\lVert{\frac{\bm{q}(\bm{p})-\tilde{\bm{q}}(\bm{p})}{\bm{q}(\bm{p})}}\right\rVert_{1},$
(10)

where $n_{r}$ is the number of discrete points in space where the PGV map
$\bm{q}$ is defined. The average mean absolute error and average mean
percentage error are obtained by average the local metrics over the discretely
sampled parameter space $\bm{P}_{x}$ and are defined via

$\displaystyle\overline{MAE}(\bm{P}_{x})$
$\displaystyle=\frac{1}{n_{x}}\sum_{i=1}^{n_{x}}{MAE}{(\bm{p}_{i})}$ (11)
$\displaystyle\overline{MAPE}(\bm{P}_{x})$
$\displaystyle=\frac{1}{n_{x}}\sum_{i=1}^{n_{x}}{MAPE}{(\bm{p}_{i})},$ (12)

where $\bm{p}_{i}\in\bm{P}_{x}$ and $n_{x}=\text{card}(\bm{P}_{x})$ is the
number of parameter instances in the dataset $\bm{P}_{x}$. When performing
cross-validation, we select the hyperparameters for each type of approximator
that result in the lowest $\overline{MAE}$ averaged across all five folds. We
then re-train the approximators on the entire training dataset using the
cross-validated hyperparameters and evaluate the $\overline{MAE}$ and
$\overline{MAPE}$, now on the completely independent testing dataset.
Altogether, this process allows us to find the best-fitting ROMs while also
giving us independent estimates of the ROM errors.

##  3 Results

###  3.1 Information content of the proper orthogonal decomposition

After generating the simulation data matrix $\bm{Q}$, we first present the
singular value distribution resulting from performing the SVD. The efficiency
and applicability of a reduced-order model for a given physical problem is
directly related to the singular value distribution. Ideally, the singular
values will rapidly decay, implying that $\bm{Q}$ has low-rank. To
characterize the singular value distribution of $\bm{Q}$ we use the relative
information content (RIC) given by

$\text{RIC}(r)=\frac{\sum_{i=1}^{r}\sigma_{i}^{2}}{\sum_{i=1}^{n_{s}}\sigma_{i}^{2}},$

where $r$ is the number of modes (left singular vectors). Using
$\text{RIC}(r)$, low-rank matrices will result in RIC $\approx 1.0$ when $r$
is small. Furthermore, the RIC may be used to guide the number of modes to use
in defining the ROM [Swischuk . (2019)], because the higher-order modes
contribute little to the solution such that they may be truncated.

In our generated datasets for all three forward models (LOH, 3D-1400, and
3D-500A) using $n_{s}=5000$, we find that only a few modes are necessary to
capture a significant portion of the information content (Figure 4).
Specifically, only 13 modes are required to recover 99% of the information
content for the data generated by LOH, 17 modes are required for 3D-1400, and
20 modes are required for 3D-500A. These results indicate that the simulated
PGV maps, independent of the forward model complexity, all have low-rank.
Since the number of modes reached to capture 99% of the information within
$\bm{Q}$ is only slightly larger for 3D velocity model setups, this suggests
that the PGV-based ROM may be only slightly more difficult to approximate
(i.e., less accurate) compared to the PGV maps associated with LOH.

Figure 4: Relative information content (RIC) plotted as a function of the
number of modes for the 5000 simulated PGV maps produced by the forward models
(LOH in red, and 3D-1400 in green, 3D-500A in blue). Note that the horizontal
axis is limited to show only the first 30 modes. We illustrate this result in
a movie (Movie S1) that shows PGV maps defined using an increasing amount of
$r$ left singular vectors, where we observe that the maps are well-
approximated by the first few modes.

###  3.2 Quantitative comparison of function approximators

Here, we present the results of training the ROMs using different types of
function approximators. In particular, we investigate the errors between the
ROM predictions and the simulated PGV maps in the testing dataset. We examine
how the ROM errors change depending on the amount of training data. We train
the ROMs with increments of $n_{s}/5$ forward models and determine the cross-
validated hyperparameters for each approximator as described in Section 2.3.
For each approximator, we show the errors (both $\overline{MAE}$ and
$\overline{MAPE}$) corresponding to the hyperparameters that result in the
lowest $\overline{MAE}$ (Figure 5). We indicate the hyperparameters that were
selected for the final evaluation of each approximator using $n_{s}$ forward
models in Table S2.

Figure 5 shows that all approximators generally result in ROMs with lower
error as more simulations are added to the training dataset. We observe that
the RBF-based ROM consistently has the lowest error out of all ROMs that we
consider, a result that is independent of the type of velocity model (either
1D or 3D). When trained with the complete set of training data, the RBF-based
ROM performs with an average absolute PGV error of $\overline{MAE}=0.07$ cm/s
for the LOH forward model, $\overline{MAE}=0.09$ cm/s for 3D-1400, and
$\overline{MAE}=0.11$ cm/s for 3D-500A. The approximate values for
$\overline{MAPE}$ are 7%, 8%, and 8%, respectively. The local slope of the
error curves, especially for $\overline{MAE}$ are variable as $n_{s}$
increases. Specifically for $n_{s}=3000$, we observe that the local slope for
$\overline{MAE}$ is approximately zero, and for all but the RBF interpolator
is actually positive. This indicates that some points in the parameter space
which are being sampled are not particularly informative for the ROM. However,
for higher values of $n_{s}$, the continuing downward slope of
$\overline{MAPE}$ for the RBF-based ROM suggests that additional training data
might further decrease the ROM errors. Similar trends are observed for all
approximators with LOH across both error metrics. Additionally, the slope of
the lines in going from $n_{s}=4000$ to 5000 is approximately the same across
all approximators. The same observations are observed for 3D-1400, with the
exception of the MLP approximator. MLP exhibits a smaller slope for large
values of $n_{s}$ in $\overline{MAE}$ and an approximately zero slope in
$\overline{MAPE}$ when $n_{s}$ increases from 4000 to 5000.

Figure 5: $\overline{MAE}$ (top row) and $\overline{MAPE}$ (bottom row) errors
for the different ROMs using the LOH (left column) and 3D-1400 (right column)
forward models. Each approximator is shown with a different colored solid line
(see legend). For each number of forward models ($n_{s}$), we evaluate the
errors using the testing dataset and the best hyperparameters selected from
5-fold cross-validation (see Table S2 for the selected hyperparameters). The
mean values of $d_{\text{nearest}}$ for the testing datasets are 0.13, 0.11,
0.10, 0.09, and 0.09 for the evaluations made using 1000, 2000, 3000, 4000,
and 5000 forward models, respectively.

In addition to examining the overall errors on the testing dataset, it is
informative to examine spatial errors within the PGV maps. We show examples of
simulated PGV maps created using all three velocity models, the predictions
from the four types of approximators, and map-views of the errors (Figure 6).
The map views of the errors illustrate locations in which the ROM produced
accurate PGV estimates, or where the PGV was over/under-predicted. For the
specific earthquake source in this example with oblique slip, the PGV is
strongest a few km northwest of the epicenter, with simulated PGVs reaching up
to 12.0 cm/s for the 3D-500A setup. Averaged across all receivers, the RBF ROM
gives the smallest error ($MAE$=0.07 cm/s for LOH), while the kNN ROM gives
the largest error ($MAE$=0.18 cm/s).

Although the errors show significant differences depending on the type of
function approximator, each approximator generally reproduces the overall
spatial distribution and amplitudes of the forward model ground motion map.
The RBF interpolator does significantly better in correctly predicting ground
motions in the areas of strongest and weakest shaking. Notably, the errors for
the RBF ROM are near zero where the shaking is the strongest. The other
approximators share some common problems, such as over-prediction at the
eastern edge of the domain and under-prediction near the southwest. The kNN,
in particular, tends to overestimate the PGVs in the near-epicentral area of
strong shaking. Overall, we conclude that the RBF ROM performs exceptionally
well at recovering the spatial distribution of ground motion compared to the
other approximators.

Figure 6: Comparison of the ROM solutions using different function
approximators for a single earthquake source (focal mechanism shown in the
first row). The second row shows the target PGV maps for LOH (left) and
3D-1400 (right). The subsequent panels show the ROM predictions and errors,
where each row is for a different approximator (see row labels at left). The
approximators are sorted from lowest to highest error (top to bottom), and
each error map is indicated with the $MAE$ in the top right corner. The error
maps outline the areas of strongest shaking where the PGV is 90% of the
largest PGV produced by each FOM.

###  3.3 ROM solutions for different sources

We next demonstrate the ROM’s ability to predict the changes in the spatial
patterns and amplitudes of peak ground motions when the hypocentral depth and
focal mechanism of the earthquake source are subject to change. We show
examples of PGV maps created with the LOH velocity model using three different
sources from the independent testing dataset (Figure 7). The PGV maps that
result from the same sources for the 3D-1400 forward model are shown in Figure
8, and the maps for 3D-500A are shown in Figure 9. In these figures, we select
three sources in $\bm{P}_{\text{test}}$ where the distance $d_{nearest}$ is
close to the expected value (0.09) of the distribution of distances for any
source in $\mathcal{P}$. We simplify the presentation by only showing the RBF
interpolator results, as this is our preferred function approximator for the
ROM. As the source focal mechanism changes, we see the ROM correctly
reproduces the changes that occur to the PGV map.

We see that the physically more complex forward models 3D-1400 and 3D-500A
produce more complex spatial PGV distributions, compared to the much smoother
distributions produced by the simpler LOH models. Comparing Figures 7 & 8, the
ROM for 3D-1400 has only slightly more difficulty in predicting the smaller-
scale features in the 3D-1400 data, and the errors are only slightly larger
compared to LOH. These small differences between LOH and 3D-1400 are
consistent with the overall trends from the testing dataset. Comparing 3D-1400
against 3D-500A, we see that 3D-500A results in larger PGVs due to the slower
wave speeds included in the velocity model. For these particular sources, the
MAE values are about 2.0 times larger for 3D-1400 compared to 3D-500A. The
complete set of error measurements on the testing dataset are shown as
histograms in Figures S4 and S5, respectively, for the three forward modeling
setups.

Figure 7: Plots of three simulated PGV maps (second column), RBF ROM
predictions (third column), and errors (fourth column) for the LOH forward
model. Each row indicates the results for a different earthquake source in the
testing dataset, where the source parameters are indicated by the focal
mechanisms in the first column. The error maps show the difference between the
FOM and ROM solutions, and indicate the mean absolute error (MAE) and distance
to the nearest training location ($d_{\text{nearest}}$). The complete set of
errors on the testing data can be viewed as histograms in Figures S4 and S5.
Figure 8: Plots of the simulated PGV maps, ROM predictions, and ROM errors for
the same sources as shown in Figure 7, but shown here for the 3D-1400 forward
model. The forward model includes 3D velocity structure (with minimum S-wave
speeds of 1400 m/s) and topography. Figure 9: Plots of the simulated PGV maps,
ROM predictions, and ROM errors for the same sources as shown in Figure 7, but
shown here for the 3D-500A forward model. The forward model includes 3D
velocity structure (with minimum S-wave speeds of 500 m/s), topography, and 3D
variable viscoelastic attenuation.

We inspect the accuracy of the ROMs for different values of the input source
parameters. In Figure 10, we show the average and standard deviation of $MAE$
for 20 evenly spaced bins for each of the four source parameters. For all
three forward models, we observe that the shallower events show larger errors
compared to deeper events. This is partially explained by the larger PGVs
expected for shallow sources, although the relative $MAPE$ errors show a non-
negligible correlation with hypocentral depth (Figure S7). There is more
variability in the errors for shallow events compared to deep events as
indicated by the standard deviations. Shallow sources introduce additional
near-source effects to ground shaking which may be underrepresented in the
training dataset following using our uniform sampling strategy. For the other
source parameters related to geometry (i.e., focal mechanism, parameterized as
strike, dip, and rake), the errors do not show significant correlations. The
standard deviations in each bin are larger for these parameters, indicating
that the ROM errors are primarily controlled by the source hypocentral depth.
Overall, these results imply that the ROM errors are less dependent on the
focal mechanism than on the source hypocentral depth.

###  3.4 Relationship between ROM errors and distance to training locations

It is reasonable to expect that the ROM errors will be smaller for the testing
parameters that have a closer neighbor in the training dataset (i.e., a
smaller value for $d_{\text{nearest}}$). We indeed see this trend for our
testing dataset (Figure S2) where the errors correlate with the distance to
the nearest neighbor. The largest error of $MAE_{i}\approx 0.8$ cm/s (for
3D-1400) occurs when the distance between the testing parameter and the
nearest training parameter is approximately equal to the fill distance of
0.17. When the Euclidean distance $d_{\text{nearest}}$ increases beyond 0.12
units, the error becomes significantly larger indicating that the
interpolation quality of the ROM degrades when the distance to the nearest
training point is too large. This relationship between $d_{\text{nearest}}$
and ROM error can be used to estimate ROM errors for earthquake sources
contained in $\mathcal{P}$ where we do not have a FOM solution. For example,
for an earthquake source with a parameter $\bm{p}$ that is 0.1 distance units
away from the nearest training FOM location, we expect the $MAE_{i}$ of the
ROM-predicted PGV map for 3D-1400 to be, on average, 0.1 cm/s.

To quantify the value added by the ROM, we examine the predicted PGV maps that
would result from a simple approach of selecting the nearest-neighboring
training FOM solution. As shown in Figure S3, the error for all forward models
is about 2.5 times larger when selecting the nearest training FOM solution,
compared to our preferred ROM estimates. This result illustrates that our ROMs
are successfully interpolating between the training FOM solutions to obtain
better estimates than could be obtained from just the simulated data.

Figure 10: Mean absolute errors for the PGV predictions on the testing
dataset, plotted as a function of each earthquake source parameter. LOH
predictions are shown in the top row, and 3D-1400predictions in the bottom
row. For each parameter, we compute the mean (black line and dots) in 20
different evenly spaced bins and plot $\pm$ 1.0 standard deviation of the
error in each bin (shaded area).

###  3.5 Quantification of the effect of depth uncertainty on PGV predictions

We now present an exemplary application using the ROMs to study ground motion
effects related to variable earthquake source parameters. In current
operational earthquake early warning algorithms, the source focal depths are
typically fixed at 8.0 or 10.0 km to minimize computational costs [Böse .
(2018), Chung . (2019)]. These depth assumptions are reasonable in California
where most earthquakes occur on shallow crustal faults [Brown . (2011)] but
may negatively affect alerting accuracy in other regions, for example,
intraslab earthquakes in the Pacific Northwest [Thompson . (2022)]. With the
ROMs, we can quantify the potential error that may be committed by an
incorrect assumption of the hypocentral depth. We use our 3D-500A ROM to make
predictions for 1 million different focal mechanisms, obtaining the PGV maps
for both a shallow (2.0 km) and a deep (8.0 km) source. The focal mechanisms
are obtained by sampling a uniform grid of 100 points along each of the three
focal mechanism parameters. Using the RBF ROM to create ground motion maps for
these 2 million earthquakes requires $\sim$1320 s on a personal
laptop111MacBook Pro (2020), Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz. This
translates to a time-to-solution for a single earthquake scenario of
$\sim$0.66 ms, or $\sim 0.6\times 10^{-6}$ CPUh222As measured by the macOS
application “Activity Monitor”.

To simplify our analysis of ground motion variability related to the source
parameters, we chose to focus on an individual site in our computational
domain. We select the Clearwater Power Plant (33.891 N, 117.609 W) near the
southeastern part of our domain (see Figure 2a). At this site, the ROM
predictions show that the PGV varies significantly depending on the focal
mechanism (Figure 11a). Averaging across all focal mechanisms, the mean
difference in PGV between the shallow and deep sources is 1.8 cm/s; however,
the maximum PGV difference is 3.6 cm/s for certain focal mechanisms. This
result implies that, in the worst-case scenario, the PGV predictions based on
an 8.0 km hypocentral depth could underestimate the true PGV at this site by
up to 3.6 cm/s if the true hypocentral depth is 2.0 km.

One advantage of the ROM approach is that we can deterministically verify
individual estimates by performing additional forward simulations. Here, we
verify the ROM-based estimate for the worst-case scenario by performing two
additional 3D-500A simulations (see Section 2.2) using the focal mechanism
that resulted in the largest ROM-predicted PGV difference (between shallow and
deep sources) at our site of interest. In the forward simulations, the PGV
difference between the shallow and deep sources is 2.8 cm/s for the site. In
comparison, the ROM-based estimate of the PGV difference is 3.6 cm/s. The
discrepancy likely results from the ROM over-estimating the PGV for the
shallow source. However, this amount of error is well within the expected
error as indicated by the mean absolute errors from the testing dataset at
this site (Figure S6). The amount of error for the ROM-estimated PGV
differences (between shallow and deep sources) that are closer to the mean
will be similar to the worst-case scenario because the errors are uncorrelated
with the focal mechanism angles (Figure 10). Overall, this experiment
highlights the potential in using ROMs to identify potentially damaging ground
motions that might be overlooked when using a smaller sample size of physics-
based simulations.

Figure 11: Histogram showing the distribution of differences between the
shallow (2.0 km) and deep (8.0 km) source PGV predictions for 1 million
different focal mechanisms (a). The ROM PGV predictions for the worst-case
focal mechanism (d) are shown for the shallow source (b) and deep source (c).
The corresponding simulations are shown in the bottom row in (e) and (f). In
the PGV maps, the site of interest is indicated by a white triangle, and the
site’s PGV value is printed.

###  3.6 Sensitivity analysis using Pearson correlation coefficients

Finally, we analyze the correlations between variable earthquake source
parameters and the PGV predictions at the same site of interest. This type of
analysis allows us to quantify the relative importance of the model inputs
(earthquake source parameters) with respect to the model output (PGV) [Homma
Saltelli (1996)]. Here, we perform a global sensitivity analysis, which is
closely related to the problem of forward uncertainty quantification
[<]e.g.,¿pagani2021enabling.

To measure the correlations, we again sample 1 million earthquakes, where the
source parameters are randomly sampled from uniform distributions. After
obtaining the ROM solutions at the Clearwater Power Plant for these sources,
we measure the correlation between each source parameter and the PGV data
using the Pearson correlation coefficient ($r_{p}$) [Pearson (1895)]. We find
a significant inverse correlation ($r_{p}=-0.65$) between the PGV predictions
and hypocentral depth (Figure 12). For the focal mechanism angles, the
correlation coefficients are near zero, implying that there is no consistent
linear relationship between any of the focal mechanism parameters and PGV.
However, there are features in the relationships between source parameters and
PGV that cannot be adequately captured by the correlation coefficient. For
example, the strike clearly affects PGV at this site by means of radiation
pattern effects which cause the visible streaking patterns. Additionally, the
kink in the PGV scaling with depth at about 5 km affects the correlation
coefficient measurement. We note that these general correlations could also be
measured from the PGV ground motions produced directly by the FOM (Figure S9),
while the 1 million ROM predictions provide finer detail compared to solely
examining the 5000 FOM measurements. Overall, while the PGV at a site of
interest may vary significantly depending on the source parameters, the
hypocentral depth is the most significant contributor to the PGV variability
in our dataset. This finding agrees with the accuracy analysis presented in
previous sections.

Lastly, we quantify our site of interest’s variability in PGV that results
from changing the source focal mechanism. We consider this type of variability
to be analogous to the inter-event terms that are often reported by GMMs which
use mixed-effects regressions, such as ASK14 [Abrahamson . (2014)]. We bin our
PGV data shown in Figure 12 into 30 equally sized hypocentral depth bins and
compute the mean and standard deviation for each bin. We find that the
standard deviation increases for increasing hypocentral depths, and ranges
from 0.18 to 0.32 natural log units (Figure S8). These variability
measurements are somewhat lower compared to GMMs, for example ASK14, which
reports an inter-event standard deviation of 0.38 natural log units for
$M_{w}=5.4$. This is an expected result, however, as the ROM-predicted PGV
variability only accounts for variations in focal mechanism and not other
source properties (e.g., stress drop) which increase the inter-event standard
deviations in GMMs. We discuss this topic further in the following section.

Figure 12: 2D density plots of the 3D-500A PGV predictions at the site of
interest, plotted against each source parameter. In each bin, the number of
earthquakes is indicated by the color bar to the right. The Pearson
correlation coefficient is indicated in the top right of each panel.

##  4 Discussion

We have constructed ROMs with the goal of predicting PGV maps for three
different model setups of increasing complexity. We indeed find that the ROMs
accurately predict PGVs (about 0.1 cm/s average error) for both the 1D layer-
over-halfspace (LOH) forward model and the 3D velocity models with topography
(3D-1400 and 3D-500A) where the ROM errors are only slightly larger for the 3D
velocity models. We find that the ROM approach works well for both linear
elastic and viscoelastic forward models of seismic wave propagation, including
low S-wave velocities, as shown by our results for the 3D-500A forward model.
The ROM ground motion errors are not uniformly distributed across the
parameter space, however, as the ROM errors tend to be larger (in both an
absolute and relative sense) for shallow events compared to deep events. While
this study uses a low-discrepancy Halton sequence to determine the FOM
training data that inform the ROM, it may be beneficial in future ROMs to non-
uniformly sample the parameter space. For example, we suspect the larger
errors for shallow sources may result from the slower near-surface velocities
which cause fewer training samples to be drawn per wavelength. As such, future
ROMs might benefit from drawing the training FOM locations proportionally to
the slowness in the velocity model, which would result in more training data
for shallow sources compared to deep sources. The error for shallow sources
could also be corrected using iterative ROM assembly [<]e.g.,¿lye2021iterative
where the errors on the testing dataset are used to guide which new parameters
will be sampled for the next iteration.

Future ROMs may investigate other function approximators beyond the four types
we considered here. We find RBF interpolation to perform the best for our
application of interest, while other studies have obtained good results using
more generalized Gaussian process regression [Guo Hesthaven (2018)]. This
approach may allow for more flexibility in the kernel length scales which can
be optimized to better fit the data, as our RBF interpolation is essentially a
special case of Gaussian process regression.

Using a Chino Hills earthquake scenario as our target simulation, we show that
ROMs can be used as on-demand ground motion simulators for modeling wave
propagation up to 1 Hz. This frequency target is consistent with the current
understanding that PGV is controlled by frequencies around 1 Hz
[<]e.g.,¿rodgersRegionalScale3D2020. However, the ROM predictions may under-
estimate real PGV for smaller earthquakes which are enriched in higher
frequency content. In addition, the minimum chosen $V_{S}$ in our forward
simulations does not account for materials with even slower wave speeds that
cause ground motion amplification in specific regions, for example, the San
Francisco Bay area muds [<]e.g.,¿hough1990sediment. To generate rapid,
physics-informed PGV maps targeted in these regions for earthquake engineering
and seismic hazard applications, the amount of computing time spent in the FOM
data collection phase would increase due to the higher cost of each forward
simulation. For example, halving the minimum shear wave velocity (or doubling
the resolved frequency) would require a factor of $2^{4}=16$ higher
computational cost (assuming spatially adaptive unstructured meshes and local
time-stepping). However, we expect the ROM approach to be able to capture the
changes to the PGV maps, as we show that there it little difference in the
rank of the simulated data matrices when we lower $V_{S,\text{min}}$ from 1400
to 500 m/s. We suspect this would also be true if the simulations resolved
higher frequencies and accounted for even slower near-surface wave speeds.
ROM-based ground motion maps from higher-frequency and physically more complex
simulations (e.g., including approximations of linear and nonlinear site
effects), may help to better understand source, path and site contributions to
ground shaking
[<]e.g.,¿hartzell1997variability,frankel2002nonlinear,hashashViscousDampingFormulation2002.

While the aim of this study is to demonstrate how the ROM approach can emulate
a ground motion simulation scenario with variable source properties, we note
in the following possible future extensions that may generalize the ROMs for
ground motion prediction of real earthquakes with arbitrary source locations
and magnitudes. For one, we expect the ROM approach will be applicable to
forward models involving more realistic earthquake source representations,
such as finite kinematic or dynamic rupture sources. In this study, we did not
parameterize forward simulations with variable moment magnitudes. While our
PGV maps could be scaled to represent different $M_{w}$ using a magnitude-
scaling term from GMMs, for example [Baltay Hanks (2014)], the point-source
approximation we use here omits finite source effects such as directivity
[Spudich . (2014)]. We note that our simulations only consider a single
epicentral location, and will not accurately predict the path effects in
regions with different velocity structures. This issue may be investigated by
extending $\mathcal{P}$ to include latitude and longitude as additional source
parameters, and we expect the ROM approach to be useful in that case to study
basin or topographic effects that depend on the source-to-site azimuth
[<]e.g.,¿[]thompsonBasinAmplificationEffects2020,stone2022topographic.
Variations to our fixed source-time function would allow us to study variable
stress drop, which has been shown to correlate with peak ground acceleration
in observational datasets [Trugman Shearer (2018)]. Recorded ground motion
observations are typically not available at the required density to build
ROMs. While we here rely on physics-based simulations to obtain PGV maps, we
expect that future work may incorporate sparse, observed ground motions into
the ROM framework to provide additional constraints on the ground motion
amplitudes and variability.

While finite earthquake source effects on PGVs may be readily included in ROMs
with additional source parameters, an alternative approach is to predict
Green’s functions instead of PGV maps. This would allow the ground motion
intensities to be computed for arbitrary source-time functions, and allow for
predictions using finite-fault sources using the superposition property of the
elastic wave equation. However, this approach might be challenging if the
Green’s functions do not possess low-rank structure like the PGV maps we
predict in this study. For this application, deep learning and neural operator
methods for time-dependent output quantities may be considered
[<]e.g.,¿[]yangSeismicWavePropagation2021a,yang2023rapid though they have yet
to be applied to 3D elastodynamic wave propagation in models that include
topography. Future directions in instantaneous ground motion estimates may
involve non-source-based approaches that use ROMs to forecast wavefield
snapshots, either used simulated or observed data [<]e.g.,¿nagata2023seismic.

##  5 Conclusion

We develop a reduced-order modeling (ROM) approach and use it to instantly
create peak ground velocity (PGV) maps that closely match computationally
expensive wave propagation simulations with variable source depth and focal
mechanisms. We show that the simulated PGV maps have low-rank structure which
motivates our use of the interpolated proper orthogonal decomposition
technique. The accuracies of the ROMs are similar when a 1D or 3D-velocity
model (including topography, viscoelastic attenuation, and slow S-wave speeds)
is used for the simulations of seismic wave propagation. We find that the ROM
errors are primarily controlled by the hypocentral depth (larger errors for
shallow sources). After comparing different approximation strategies with
cross-validation, we find that radial basis function interpolation performs
well, showing small (mean absolute PGV error $<0.1$ cm/s) errors on an
independent testing dataset. We lastly demonstrate the potential of using ROMs
to answer questions about ground motion variability that would ordinarily be
computationally intractable when using full complexity forward simulations of
wave propagation. We conclude that the reduced-order modeling approach defined
via the interpolated proper orthogonal decomposition is a viable technique to
obtain physics-informed estimates of the relationship between source
properties and PGVs for ground motion simulation scenarios.

##  6 Open Research

###  6.1 Data Availability Statement

We use the Scikit-learn [Pedregosa . (2011)] and SciPy [Virtanen . (2020)]
Python modules for their implementations of the approximators used in this
study. We created the figures for this work using matplotlib [Hunter (2007)],
cartopy [Met Office (2010)], and obspy [Beyreuther . (2010)]. We provide the
simulated PGV maps that we used for training the ROMs at the public GitHub
repository for this paper [Rekoske . (2023)]. This repository provides a
Jupyter notebook which demonstrates using the iPOD reduced-order modeling
technique to create ROMs for predicting new PGV maps.

###### Acknowledgements.

We would like to thank the Editor Rachel Abercrombie and the Associate Editor
for their comments that helped to significantly improve this manuscript. We
also greatly appreciate a review from Jack Muir who helped us improve the
quality of the forward simulations and clarify our methodology, as well as an
anonymous reviewer whose feedback helped us improve our presentation of the
ROM error analysis and clarify the benefits of our reduced-order model. This
material is based upon work supported by the National Science Foundation
Graduate Research Fellowship Program under Grant No. DGE-2038238. Any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect the views
of the National Science Foundation. We gratefully acknowledge the Gauss Centre
for Supercomputing e.V. (www.gauss-centre.eu) for providing computing time on
the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre
(www.lrz.de), in project pn49ha. We acknowledge additional support from the
National Science Foundation (NSF Grant No. EAR-2121568), the National
Aeronautics and Space Administration (80NSSC20K0495), from the European
Union’s Horizon 2020 Research and Innovation Programme (TEAR ERC Starting,
grant No. 852992), Horizon Europe (DT-GEO, grant No. 101058129, Geo-INQUIRE,
grant No.101058518, and ChEESE-2P, grant No. 101093038), and the Southern
California Earthquake Center (SCEC award 22135).

## References

  * Aagaard . (2001) aagaardCharacterizationNearsourceGround2001Aagaard, BT., Hall, JF. Heaton, TH.  2001\.  Characterization of Near-Source Ground Motions with Earthquake Simulations Characterization of near-source ground motions with earthquake simulations.  Earthquake Spectra172177–207.
  * Abrahamson . (2014) abrahamson2014summaryAbrahamson, NA., Silva, WJ. Kamai, R.  2014\.  Summary of the ASK14 ground motion relation for active crustal regions Summary of the ASK14 ground motion relation for active crustal regions.  Earthquake Spectra3031025–1055.
  * Anderson Brune (1999) andersonProbabilisticSeismicHazard1999Anderson, JG. Brune, JN.  1999\.  Probabilistic Seismic Hazard Analysis without the Ergodic Assumption Probabilistic seismic hazard analysis without the ergodic assumption.  Seismological Research Letters70119–28.
  * Audouze . (2009) audouzeReducedorderModelingParameterized2009Audouze, C., De Vuyst, F. Nair, P.  2009\.  Reduced-Order Modeling of Parameterized PDEs Using Time–Space-Parameter Principal Component Analysis Reduced-order modeling of parameterized PDEs using time–space-parameter principal component analysis.  International journal for numerical methods in engineering8081025–1057.
  * Bailey . (2022) bailey2022adaptingBailey, MD., Bandyopadhyay, S. Nychka, D.  2022\.  Adapting conditional simulation using circulant embedding for irregularly spaced spatial data Adapting conditional simulation using circulant embedding for irregularly spaced spatial data.  Stat111e446.
  * Baker . (2014) Baker2014Baker, JW., Luco, N., Abrahamson, NA., Graves, RW., Maechling, PJ. Olsen, KB.  2014\.  Engineering uses of physics-based ground motion simulations Engineering uses of physics-based ground motion simulations.  Proceedings of the Tenth US Conference on Earthquake Engineering. Proceedings of the tenth us conference on earthquake engineering.
  * Baltay Hanks (2014) baltay2014understandingBaltay, AS. Hanks, TC.  2014\.  Understanding the magnitude dependence of PGA and PGV in NGA-West 2 data Understanding the magnitude dependence of pga and pgv in nga-west 2 data.  Bulletin of the Seismological Society of America10462851–2865.
  * Baltay . (2017) baltayUncertaintyVariabilityEarthquake2017Baltay, AS., Hanks, TC. Abrahamson, NA.  2017\.  Uncertainty, Variability, and Earthquake Physics in Ground-Motion Prediction Equations Uncertainty, variability, and earthquake physics in ground-motion prediction equations.  Bulletin of the Seismological Society of America10741754–1772.
  * Benner . (2015) bennerSurveyProjectionbasedModel2015Benner, P., Gugercin, S. Willcox, K.  2015\.  A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM review574483–531.
  * Berkooz . (1993) berkoozProperOrthogonalDecomposition1993Berkooz, G., Holmes, P. Lumley, JL.  1993\.  The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows The proper orthogonal decomposition in the analysis of turbulent flows.  Annual review of fluid mechanics251539–575.
  * Beyreuther . (2010) beyreutherObsPyPythonToolbox2010Beyreuther, M., Barsch, R., Krischer, L., Megies, T., Behr, Y. Wassermann, J.  2010\.  ObsPy: A Python Toolbox for Seismology ObsPy: A Python toolbox for seismology.  Seismological Research Letters813530–533.
  * Bommer . (2010) bommerSelectionGroundmotionPrediction2010Bommer, JJ., Douglas, J., Scherbaum, F., Cotton, F., Bungum, H. Fah, D.  2010\.  On the Selection of Ground-Motion Prediction Equations for Seismic Hazard Analysis On the selection of ground-motion prediction equations for seismic hazard analysis.  Seismological Research Letters815783–793.
  * Boore (1983) booreStochasticSimulationHighfrequency1983Boore, DM.  1983\.  Stochastic Simulation of High-Frequency Ground Motions Based on Seismological Models of the Radiated Spectra Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra.  Bulletin of the Seismological Society of America736A1865–1894.  10.1785/BSSA07306A1865
  * Boore (2010) booreOrientationindependentNongeometricmeanMeasures2010Boore, DM.  2010\.  Orientation-Independent, Nongeometric-Mean Measures of Seismic Intensity from Two Horizontal Components of Motion Orientation-independent, nongeometric-mean measures of seismic intensity from two horizontal components of motion.  Bulletin of the Seismological Society of America10041830–1835.  10.1785/0120090400
  * Böse . (2018) boseFinDerImprovedRealtime2018bBöse, M., Smith, DE., Felizardo, C., Meier, MA., Heaton, TH. Clinton, JF.  2018\.  FinDer v.2: Improved Real-Time Ground-Motion Predictions for M2–M9 with Seismic Finite-Source Characterization FinDer v.2: Improved real-time ground-motion predictions for M2–M9 with seismic finite-source characterization.  Geophysical Journal International2121725–742.  10.1093/gji/ggx430
  * Bozorgnia . (2021) bozorgniaNGASubductionResearchProgram2021Bozorgnia, Y., Abrahamson, NA., Ahdi, SK., Ancheta, TD., Atik, LA., Archuleta, RJ.SJ Chiou, B.  2021\.  NGA-Subduction Research Program NGA-Subduction research program.  Earthquake Spectra382783–798.
  * Bozorgnia . (2014) bozorgniaNGAWest2ResearchProject2014Bozorgnia, Y., Abrahamson, NA., Atik, LA., Ancheta, TD., Atkinson, GM., Baker, JW.Chiou, BSJ.  2014\.  NGA-West2 Research Project NGA-West2 research project.  Earthquake Spectra303973–987.
  * Breiman (2001) breimanRandomForests2001Breiman, L.  2001\.  Random Forests Random forests.  Machine learning4515–32.
  * Brown . (2011) brown2011developmentBrown, HM., Allen, RM., Hellweg, M., Khainovski, O., Neuhauser, D. Souf, A.  2011\.  Development of the ElarmS methodology for earthquake early warning: Realtime application in California and offline testing in Japan Development of the elarms methodology for earthquake early warning: Realtime application in California and offline testing in Japan.  Soil Dynamics and Earthquake Engineering312188–200.
  * Brune (1970) bruneTectonicStressSpectra1970Brune, JN.  1970\.  Tectonic Stress and the Spectra of Seismic Shear Waves from Earthquakes Tectonic stress and the spectra of seismic shear waves from earthquakes.  Journal of Geophysical Research75264997–5009.
  * Bui-Thanh . (2003) bui2003properBui-Thanh, T., Damodaran, M. Willcox, K.  2003\.  Proper orthogonal decomposition extensions for parametric applications in compressible aerodynamics Proper orthogonal decomposition extensions for parametric applications in compressible aerodynamics.  21st AIAA applied aerodynamics conference 21st aiaa applied aerodynamics conference ( 4213).
  * Chung . (2019) chungOptimizingEarthquakeEarly2019Chung, AI., Henson, I. Allen, RM.  2019\.  Optimizing Earthquake Early Warning Performance: ElarmS-3 Optimizing earthquake early warning performance: ElarmS-3.  Seismological Research Letters902A727–743.  10.1785/0220180192
  * Cochran . (2019) cochran2019eventCochran, ES., Bunn, J., Minson, SE., Baltay, AS., Kilb, DL., Kodera, Y. Hoshiba, M.  2019\.  Event detection performance of the PLUM earthquake early warning algorithm in southern California Event detection performance of the PLUM earthquake early warning algorithm in southern California.  Bulletin of the Seismological Society of America10941524–1541.
  * Cui . (2013) Cui2013Cui, Y., Poyraz, E., Olsen, KB., Zhou, J., Withers, K., Callaghan, S.others  2013\.  Physics-based seismic hazard analysis on petascale heterogeneous supercomputers Physics-based seismic hazard analysis on petascale heterogeneous supercomputers.  SC’13: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis Sc’13: Proceedings of the international conference on high performance computing, networking, storage and analysis ( 1–12).
  * De Marchi Schaback (2010) de2010stabilityDe Marchi, S. Schaback, R.  2010\.  Stability of kernel-based interpolation Stability of kernel-based interpolation.
  * Druault . (2005) druaultUseProperOrthogonal2005Druault, P., Guibert, P. Alizon, F.  2005\.  Use of Proper Orthogonal Decomposition for Time Interpolation from PIV Data Use of proper orthogonal decomposition for time interpolation from PIV data.  Experiments in Fluids3961009–1023.  10.1007/s00348-005-0035-3
  * Dumbser Käser (2006) dumbserArbitraryHighorderDiscontinuous2006Dumbser, M. Käser, M.  2006\.  An Arbitrary High-Order Discontinuous Galerkin Method for Elastic Waves on Unstructured Meshes – II. The Three-Dimensional Isotropic Case An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes – II. The three-dimensional isotropic case.  Geophysical Journal International1671319–336.  10.1111/j.1365-246X.2006.03120.x
  * Engler . (2022) engler2022partitioningEngler, DT., Worden, CB., Thompson, EM. Jaiswal, KS.  2022\.  Partitioning Ground Motion Uncertainty When Conditioned on Station Data Partitioning ground motion uncertainty when conditioned on station data.  Bulletin of the Seismological Society of America11221060–1079.
  * Field . (2014) field2014uniformField, EH., Arrowsmith, RJ., Biasi, GP., Bird, P., Dawson, TE., Felzer, KR.others  2014\.  Uniform California earthquake rupture forecast, version 3 (UCERF3)—The time-independent model Uniform california earthquake rupture forecast, version 3 (ucerf3)—the time-independent model.  Bulletin of the Seismological Society of America10431122–1180.
  * A. Frankel . (2018) frankelBroadbandSyntheticSeismograms2018Frankel, A., Wirth, E., Marafi, N., Vidale, J. Stephenson, W.  2018\.  Broadband Synthetic Seismograms for Magnitude 9 Earthquakes on the Cascadia Megathrust Nased on 3D Simulations and Stochastic Synthetics, Part 1: Methodology and Overall Results Broadband synthetic seismograms for magnitude 9 earthquakes on the Cascadia megathrust nased on 3D simulations and stochastic synthetics, part 1: Methodology and overall results.  Bulletin of the Seismological Society of America1085A2347–2369.
  * AD. Frankel . (2002) frankel2002nonlinearFrankel, AD., Carver, DL. Williams, RA.  2002\.  Nonlinear and linear site response and basin effects in Seattle for the M 6.8 Nisqually, Washington, earthquake Nonlinear and linear site response and basin effects in seattle for the m 6.8 nisqually, washington, earthquake.  Bulletin of the Seismological Society of America9262090–2109.
  * Furumura . (2019) furumura2019earlyFurumura, T., Maeda, T. Oba, A.  2019\.  Early forecast of long-period ground motions via data assimilation of observed ground motions and wave propagation simulations Early forecast of long-period ground motions via data assimilation of observed ground motions and wave propagation simulations.  Geophysical Research Letters461138–147.
  * Gallovič . (2019) gallovicBayesianDynamicFinitefault2019bGallovič, F., Valentová, Ľ., Ampuero, JP. Gabriel, AA.  2019\.  Bayesian Dynamic Finite-Fault Inversion: 1. Method and Synthetic Test Bayesian dynamic finite-fault inversion: 1. method and synthetic test.  Journal of Geophysical Research: Solid Earth12476949–6969.  10.1029/2019JB017510
  * Geuzaine Remacle (2009) geuzaineGmsh3DFinite2009Geuzaine, C. Remacle, JF.  2009\.  Gmsh: A 3-D Finite Element Mesh Generator with Built-in Pre-and Post-Processing Facilities Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities.  International journal for numerical methods in engineering79111309–1331.
  * Goulet . (2021) gouletNGAEastGroundMotionCharacterization2021Goulet, CA., Bozorgnia, Y., Kuehn, N., Al Atik, L., Youngs, RR., Graves, RW. Atkinson, GM.  2021\.  NGA-East Ground-Motion Characterization Model Part I: Summary of Products and Model Development NGA-East ground-motion characterization model part I: Summary of products and model development.  Earthquake Spectra371_suppl1231–1282.  10.1177/87552930211018723
  * Graves . (2011) gravesCyberShakePhysicsbasedSeismic2011Graves, R., Jordan, TH., Callaghan, S., Deelman, E., Field, E., Juve, G.Vahi, K.  2011\.  CyberShake: A Physics-Based Seismic Hazard Model for Southern California CyberShake: A physics-based seismic hazard model for Southern California.  Pure and Applied Geophysics1683367–381.  10.1007/s00024-010-0161-6
  * Graves Pitarka (2016) gravesKinematicGroundmotionSimulations2016Graves, R. Pitarka, A.  2016\.  Kinematic Ground-Motion Simulations on Rough Faults Including Effects of 3D Stochastic Velocity Perturbations Kinematic ground-motion simulations on rough faults including effects of 3D stochastic velocity perturbations.  Bulletin of the Seismological Society of America10652136–2153.
  * Guo Hesthaven (2018) guoReducedOrderModeling2018Guo, M. Hesthaven, JS.  2018\.  Reduced Order Modeling for Nonlinear Structural Analysis Using Gaussian Process Regression Reduced order modeling for nonlinear structural analysis using Gaussian process regression.  Computer methods in applied mechanics and engineering341807–826.
  * Halton (1960) haltonEfficiencyCertainQuasirandom1960Halton, JH.  1960\.  On the Efficiency of Certain Quasi-Random Sequences of Points in Evaluating Multi-Dimensional Integrals On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals.  Numerische Mathematik2184–90.
  * Hartzell . (1997) hartzell1997variabilityHartzell, S., Cranswick, E., Frankel, A., Carver, D. Meremonte, M.  1997\.  Variability of site response in the Los Angeles urban area Variability of site response in the los angeles urban area.  Bulletin of the Seismological Society of America8761377–1400.
  * Hashash Park (2002) hashashViscousDampingFormulation2002Hashash, YM. Park, D.  2002\.  Viscous Damping Formulation and High Frequency Motion Propagation in Non-Linear Site Response Analysis Viscous damping formulation and high frequency motion propagation in non-linear site response analysis.  Soil Dynamics and Earthquake Engineering227611–624.
  * Heinecke . (2014) Heinecke2014Heinecke, A., Breuer, A., Rettenberger, S., Bader, M., Gabriel, AA., Pelties, C.others  2014\.  Petascale high order dynamic rupture earthquake simulations on heterogeneous supercomputers Petascale high order dynamic rupture earthquake simulations on heterogeneous supercomputers.  SC’14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis Sc’14: Proceedings of the international conference for high performance computing, networking, storage and analysis ( 3–14).
  * Hesthaven . (2022) hesthavenReducedBasisMethods2022Hesthaven, JS., Pagliantini, C. Rozza, G.  2022\.  Reduced Basis Methods for Time-Dependent Problems Reduced basis methods for time-dependent problems.  Acta Numerica31265–345.
  * Hesthaven Ubbiali (2018) hesthavenNonintrusiveReducedOrder2018bHesthaven, JS. Ubbiali, S.  2018\.  Non-Intrusive Reduced Order Modeling of Nonlinear Problems Using Neural Networks Non-intrusive reduced order modeling of nonlinear problems using neural networks.  Journal of Computational Physics36355–78.  10.1016/j.jcp.2018.02.037
  * Hinton (1990) hintonConnectionistLearningProcedures1990Hinton, GE.  1990\.  Connectionist Learning Procedures Connectionist learning procedures.  Machine Learning Machine learning ( 555–610).  Elsevier.
  * Homma Saltelli (1996) homma1996importanceHomma, T. Saltelli, A.  1996\.  Importance measures in global sensitivity analysis of nonlinear models Importance measures in global sensitivity analysis of nonlinear models.  Reliability Engineering & System Safety5211–17.
  * Hoshiba Aoki (2015) hoshiba2015numericalHoshiba, M. Aoki, S.  2015\.  Numerical Shake Prediction for Earthquake Early Warning: Data Assimilation, Real-Time Shake Mapping, and Simulation of Wave Propagation Numerical shake prediction for earthquake early warning: Data assimilation, real-time shake mapping, and simulation of wave propagation.  Bulletin of the Seismological Society of America10531324–1338.
  * Hough . (1990) hough1990sedimentHough, SE., Friberg, P., Busby, R., Field, E., Jacob, K. Borcherdt, R.  1990\.  Sediment-induced amplification and the collapse of the Nimitz Freeway Sediment-induced amplification and the collapse of the Nimitz Freeway.  Nature3446269853–855.
  * Hu . (2022) huHzDeterministic3D2022Hu, Z., Olsen, KB. Day, SM.  2022\.  0–5 Hz Deterministic 3-D Ground Motion Simulations for the 2014 La Habra, California, Earthquake 0–5 Hz deterministic 3-D ground motion simulations for the 2014 La Habra, California, Earthquake.  Geophysical Journal International23032162–2182.
  * Hunter (2007) hunterMatplotlib2DGraphics2007Hunter, JD.  2007\.  Matplotlib: A 2D Graphics Environment Matplotlib: A 2D graphics environment.  Computing in science & engineering90390–95.
  * Julian . (1998) julian1998nonJulian, BR., Miller, AD. Foulger, G.  1998\.  Non-double-couple earthquakes 1. Theory Non-double-couple earthquakes 1. theory.  Reviews of Geophysics364525–549.
  * Käser . (2008) Kaeser2008Käser, M., Hermann, V. Puente, Jdl.  2008\.  Quantitative accuracy analysis of the discontinuous Galerkin method for seismic wave propagation Quantitative accuracy analysis of the discontinuous galerkin method for seismic wave propagation.  Geophysical Journal International1733990–999.
  * Kingma Ba (2017) kingmaAdamMethodStochastic2017Kingma, DP. Ba, J.  2017\.  Adam: A Method for Stochastic Optimization Adam: A Method for Stochastic Optimization ( arXiv:1412.6980).  arXiv.
  * Kodera . (2018) kodera2018propagationKodera, Y., Yamada, Y., Hirano, K., Tamaribuchi, K., Adachi, S., Hayashimoto, N.Hoshiba, M.  2018\.  The Propagation of Local Undamped Motion (PLUM) Method: A Simple and Robust Seismic Wavefield Estimation Approach for Earthquake Early Warning The propagation of local undamped motion (PLUM) method: A simple and robust seismic wavefield estimation approach for earthquake early warning.  Bulletin of the Seismological Society of America1082983–1003.
  * Kohler . (2017) kohlerEarthquakeEarlyWarning2017Kohler, MD., Cochran, ES., Given, D., Guiwits, S., Neuhauser, D., Henson, I.Schwarz, S.  2017\.  Earthquake Early Warning ShakeAlert System: West Coast Wide Production Prototype Earthquake early warning ShakeAlert system: West coast wide production prototype.  Seismological Research Letters89199–107.  10.1785/0220170140
  * Komatitsch Tromp (1999) komatitschIntroductionSpectralElement1999aKomatitsch, D. Tromp, J.  1999\.  Introduction to the Spectral Element Method for Three-Dimensional Seismic Wave Propagation Introduction to the spectral element method for three-dimensional seismic wave propagation.  Geophysical journal international1393806–822.
  * Krenz . (2021) Krenz2021Krenz, L., Uphoff, C., Ulrich, T., Gabriel, AA., Abrahams, L., Dunham, E. Bader, M.  2021\.  3D acoustic-elastic coupling with gravity: the dynamics of the 2018 Palu, Sulawesi earthquake and tsunami 3D acoustic-elastic coupling with gravity: the dynamics of the 2018 Palu, Sulawesi earthquake and tsunami.  Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. Proceedings of the international conference for high performance computing, networking, storage and analysis.  New York, NY, USAACM.  https://doi.org/10.1145/3458817.3476173
  * Lazzaro Montefusco (2002) lazzaroRadialBasisFunctions2002Lazzaro, D. Montefusco, LB.  2002\.  Radial Basis Functions for the Multivariate Interpolation of Large Scattered Data Sets Radial basis functions for the multivariate interpolation of large scattered data sets.  Journal of Computational and Applied Mathematics1401-2521–536.
  * Lye . (2021) lye2021iterativeLye, KO., Mishra, S., Ray, D. Chandrashekar, P.  2021\.  Iterative surrogate model optimization (ISMO): an active learning algorithm for PDE constrained optimization with deep neural networks Iterative surrogate model optimization (ismo): an active learning algorithm for pde constrained optimization with deep neural networks.  Computer Methods in Applied Mechanics and Engineering374113575.
  * Meletti . (2021) melettiNewItalianSeismic2021Meletti, C., Marzocchi, W., D’amico, V., Lanzano, G., Luzi, L., Martinelli, F.Visini, F.  2021\.  The New Italian Seismic Hazard Model (MPS19) The new Italian seismic hazard model (MPS19).  Annals of Geophysics641.
  * Met Office (2010) metofficeCartopyCartographicPython2010Met Office, UK.  2010\.  Cartopy: A Cartographic Python Library with a Matplotlib Interface Cartopy: A cartographic python library with a matplotlib interface.  Exeter, Devon.
  * Milner . (2021) milnerPhysicsbasedNonergodicPSHA2021Milner, KR., Shaw, BE., Goulet, CA., Richards-Dinger, KB., Callaghan, S., Jordan, TH.Field, EH.  2021\.  Toward Physics-Based Nonergodic PSHA: A Prototype Fully Deterministic Seismic Hazard Model for Southern California Toward physics-based nonergodic PSHA: A prototype fully deterministic seismic hazard model for southern California.  Bulletin of the Seismological Society of America1112898–915.
  * Moschetti . (2017) moschetti3DGroundmotionSimulations2017Moschetti, MP., Hartzell, S., Ramírez-Guzmán, L., Frankel, AD., Angster, SJ. Stephenson, WJ.  2017\.  3D Ground-Motion Simulations of Mw 7 Earthquakes on the Salt Lake City Segment of the Wasatch Fault Zone: Variability of Long-Period (T $\geq$ 1 s) Ground Motions and Sensitivity to Kinematic Rupture Parameters 3D ground-motion simulations of Mw 7 earthquakes on the Salt Lake City segment of the Wasatch fault zone: Variability of long-period (T $\geq$ 1 s) ground motions and sensitivity to kinematic rupture parameters.  Bulletin of the Seismological Society of America10741704–1723.
  * Moseley . (2018) moseleyFastApproximateSimulation2018Moseley, B., Markham, A. Nissen-Meyer, T.  2018\.  Fast Approximate Simulation of Seismic Waves with Deep Learning Fast approximate simulation of seismic waves with deep learning.  arXiv preprint arXiv:1807.06873.
  * Moseley . (2020) moseleyDeepLearningFast2020Moseley, B., Nissen-Meyer, T. Markham, A.  2020\.  Deep Learning for Fast Simulation of Seismic Waves in Complex Media Deep learning for fast simulation of seismic waves in complex media.  Solid Earth1141527–1549.  10.5194/se-11-1527-2020
  * Nagata . (2023) nagata2023seismicNagata, T., Nakai, K., Yamada, K., Saito, Y., Nonomura, T., Kano, M.Nagao, H.  2023\.  Seismic wavefield reconstruction based on compressed sensing using data-driven reduced-order model Seismic wavefield reconstruction based on compressed sensing using data-driven reduced-order model.  Geophysical Journal International233133–50.
  * Oba . (2020) oba2020dataOba, A., Furumura, T. Maeda, T.  2020\.  Data assimilation-based early forecasting of long-period ground motions for large earthquakes along the Nankai trough Data assimilation-based early forecasting of long-period ground motions for large earthquakes along the Nankai trough.  Journal of Geophysical Research: Solid Earth1256e2019JB019047.
  * K. Olsen . (2003) olsen2003estimationOlsen, K., Day, S. Bradley, C.  2003\.  Estimation of Q for long-period (¿ 2 sec) waves in the Los Angeles basin Estimation of q for long-period (¿ 2 sec) waves in the los angeles basin.  Bulletin of the Seismological Society of America932627–638.
  * KB. Olsen . (2006) olsenStrongShakingAngeles2006Olsen, KB., Day, SM., Minster, JB., Cui, Y., Chourasia, A., Faerman, M.Jordan, T.  2006\.  Strong Shaking in Los Angeles Expected from Southern San Andreas Earthquake Strong shaking in Los Angeles expected from southern San Andreas earthquake.  Geophysical Research Letters337.  10.1029/2005GL025472
  * Pagani Manzoni (2021) pagani2021enablingPagani, S. Manzoni, A.  2021\.  Enabling forward uncertainty quantification and sensitivity analysis in cardiac electrophysiology by reduced order modeling and machine learning Enabling forward uncertainty quantification and sensitivity analysis in cardiac electrophysiology by reduced order modeling and machine learning.  International Journal for Numerical Methods in Biomedical Engineering376e3450.
  * Parker Baltay (2022) parkerEmpiricalMapBasedNonergodic2022Parker, GA. Baltay, AS.  2022\.  Empirical Map-Based Nonergodic Models of Site Response in the Greater Los Angeles Area Empirical map-based nonergodic models of site response in the Greater Los Angeles area.  Bulletin of the Seismological Society of America11231607–1629.
  * Pearson (1895) pearson1895viiPearson, K.  1895\.  VII. Note on regression and inheritance in the case of two parents Vii. note on regression and inheritance in the case of two parents.  proceedings of the royal society of London58347-352240–242.
  * Pedregosa . (2011) pedregosaScikitlearnMachineLearning2011Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O.Duchesnay, É.  2011\.  Scikit-Learn: Machine Learning in Python Scikit-learn: Machine Learning in Python.  Journal of Machine Learning Research12852825–2830.
  * Petersen . (2020) petersen2018UpdateUS2020Petersen, MD., Shumway, AM., Powers, PM., Mueller, CS., Moschetti, MP., Frankel, AD.Boyd, OS.  2020\.  The 2018 Update of the US National Seismic Hazard Model: Overview of Model and Implications The 2018 update of the US National Seismic Hazard Model: Overview of model and implications.  Earthquake Spectra3615–41.
  * Pitarka . (2022) pitarkaDeterministic3DGroundMotion2022bPitarka, A., Akinci, A., De Gori, P. Buttinelli, M.  2022\.  Deterministic 3D Ground-Motion Simulations (0–5 Hz) and Surface Topography Effects of the 30 October 2016 M w 6.5 Norcia, Italy, Earthquake Deterministic 3D Ground-Motion Simulations (0–5 Hz) and Surface Topography Effects of the 30 October 2016 M w 6.5 Norcia, Italy, Earthquake.  Bulletin of the Seismological Society of America1121262–286.
  * Rasht-Behesht . (2022) rasht-beheshtPhysicsInformedNeuralNetworks2022Rasht-Behesht, M., Huber, C., Shukla, K. Karniadakis, GE.  2022\.  Physics-Informed Neural Networks (PINNs) for Wave Propagation and Full Waveform Inversions Physics-informed neural networks (PINNs) for wave propagation and full waveform inversions.  Journal of Geophysical Research: Solid Earth1275e2021JB023120.
  * Rekoske . (2023) rekoske_john_2023_8118300Rekoske, J., Gabriel, AA. May, D.  2023\.  jrekoske/Instantaneous-PhysicsBased- GroundMotionMaps: Publication release (Version 1.0.0) jrekoske/Instantaneous-PhysicsBased- GroundMotionMaps: Publication release (Version 1.0.0) [Software].  Zenodo.  https://doi.org/10.5281/zenodo.8170242 10.5281/zenodo.8170242
  * Rettenberger . (2016) Rettenberger2016Rettenberger, S., Meister, O., Bader, M. Gabriel, AA.  2016\.  ASAGI: A Parallel Server for Adaptive Geoinformation Asagi: A parallel server for adaptive geoinformation.  Proceedings of the Exascale applications and Software Conference 2016 Proceedings of the exascale applications and software conference 2016 ( 2:1-2:9).  New York, NY, USAACM.  https://doi.org/10.1145/2938615.2938618
  * Ripperger . (2008) rippergerVariabilityNearfieldGround2008Ripperger, J., Mai, PM. Ampuero, JP.  2008\.  Variability of Near-Field Ground Motion from Dynamic Earthquake Rupture Simulations Variability of near-field ground motion from dynamic earthquake rupture simulations.  Bulletin of the Seismological Society of America9831207–1228.
  * Rodgers . (2020) rodgersRegionalScale3D2020Rodgers, AJ., Pitarka, A., Pankajakshan, R., Sjögreen, B. Petersson, NA.  202008\.  Regional-Scale 3D Ground-Motion Simulations of Mw 7 Earthquakes on the Hayward Fault, Northern California Resolving Frequencies 0–10 Hz and Including Site-Response Corrections Regional-Scale 3D Ground-Motion Simulations of Mw 7 Earthquakes on the Hayward Fault, Northern California Resolving Frequencies 0–10 Hz and Including Site-Response Corrections.  Bulletin of the Seismological Society of America11062862–2881.  10.1785/0120200147
  * Sahakian . (2019) sahakianGroundMotionResiduals2019Sahakian, VJ., Baltay, A., Hanks, TC., Buehler, J., Vernon, FL., Kilb, D. Abrahamson, NA.  2019\.  Ground Motion Residuals, Path Effects, and Crustal Properties: A Pilot Study in Southern California Ground motion residuals, path effects, and crustal properties: A pilot study in Southern California.  Journal of Geophysical Research: Solid Earth12465738–5753.  10.1029/2018JB016796
  * Savran Olsen (2020) Savran2020Savran, W. Olsen, K.  2020\.  Kinematic rupture generator based on 3-D spontaneous rupture simulations along geometrically rough faults Kinematic rupture generator based on 3-d spontaneous rupture simulations along geometrically rough faults.  Journal of Geophysical Research: Solid Earth12510e2020JB019464.
  * SCEDC (2013) SCEDC2013SCEDC.  2013\.  Southern California Earthquake Data Center. Southern California Earthquake Data Center.  https://scedc.caltech.edu https://doi.org/10.7909/C3WD3xH1
  * Shaw . (2015) shawUnifiedStructuralRepresentation2015Shaw, JH., Plesch, A., Tape, C., Suess, MP., Jordan, TH., Ely, G.Munster, J.  2015\.  Unified Structural Representation of the Southern California Crust and Upper Mantle Unified structural representation of the Southern California crust and upper mantle.  Earth and Planetary Science Letters4151–15.  10.1016/j.epsl.2015.01.016
  * Small . (2017) smallSCECUnifiedCommunity2017aSmall, P., Gill, D., Maechling, PJ., Taborda, R., Callaghan, S., Jordan, TH.Goulet, C.  2017\.  The SCEC Unified Community Velocity Model Software Framework The SCEC unified community velocity model software framework.  Seismological Research Letters8861539–1552.
  * Smith . (2020) smithEikonetSolvingEikonal2020Smith, JD., Azizzadenesheli, K. Ross, ZE.  2020\.  Eikonet: Solving the Eikonal Equation with Deep Neural Networks Eikonet: Solving the eikonal equation with deep neural networks.  IEEE Transactions on Geoscience and Remote Sensing591210685–10696.
  * Spudich . (2014) spudich2014comparisonSpudich, P., Rowshandel, B., Shahi, SK., Baker, JW. Chiou, BSJ.  2014\.  Comparison of NGA-West2 directivity models Comparison of NGA-West2 directivity models.  Earthquake Spectra3031199–1221.
  * Stirling . (2012) stirlingNationalSeismicHazard2012Stirling, M., McVerry, G., Gerstenberger, M., Litchfield, N., Van Dissen, R., Berryman, K.Jacobs, K.  2012\.  National Seismic Hazard Model for New Zealand: 2010 Update National Seismic Hazard Model for New Zealand: 2010 Update.  Bulletin of the Seismological Society of America10241514–1542.  10.1785/0120110170
  * Stone . (2022) stone2022topographicStone, I., Wirth, EA. Frankel, AD.  2022\.  Topographic Response to Simulated Mw 6.5–7.0 Earthquakes on the Seattle Fault Topographic response to simulated Mw 6.5–7.0 earthquakes on the Seattle Fault.  Bulletin of the Seismological Society of America11231436–1462.
  * Swischuk . (2019) swischukProjectionbasedModelReduction2019Swischuk, R., Mainini, L., Peherstorfer, B. Willcox, K.  2019\.  Projection-Based Model Reduction: Formulations for Physics-Based Machine Learning Projection-based model reduction: Formulations for physics-based machine learning.  Computers & Fluids179704–717.  10.1016/j.compfluid.2018.07.021
  * Taborda Bielak (2013) tabordaGroundmotionSimulationValidation2013Taborda, R. Bielak, J.  2013\.  Ground-Motion Simulation and Validation of the 2008 Chino Hills, California, Earthquake Ground-motion simulation and validation of the 2008 Chino Hills, California, earthquake.  Bulletin of the Seismological Society of America1031131–156.
  * Tamhidi . (2021) tamhidiConditionedSimulationGround2021Tamhidi, A., Kuehn, N., Ghahari, SF., Rodgers, AJ., Kohler, MD., Taciroglu, E. Bozorgnia, Y.  202109\.  Conditioned Simulation of Ground-Motion Time Series at Uninstrumented Sites Using Gaussian Process Regression Conditioned Simulation of Ground-Motion Time Series at Uninstrumented Sites Using Gaussian Process Regression.  Bulletin of the Seismological Society of America1121331–347.  10.1785/0120210054
  * Taufiqurrahman . (2022) taufiqurrahmanBroadbandDynamicRupture2022aTaufiqurrahman, T., Gabriel, AA., Ulrich, T., Valentová, L. Gallovič, F.  2022\.  Broadband Dynamic Rupture Modeling With Fractal Fault Roughness, Frictional Heterogeneity, Viscoelasticity and Topography: The 2016 Mw 6.2 Amatrice, Italy Earthquake Broadband Dynamic Rupture Modeling With Fractal Fault Roughness, Frictional Heterogeneity, Viscoelasticity and Topography: The 2016 Mw 6.2 Amatrice, Italy Earthquake.  Geophysical Research Letters4922e2022GL098872.  10.1029/2022GL098872
  * Thompson . (2022) thompsonEffectFixingEarthquake2022Thompson, M., Hartog, JR. Wirth, EA.  2022\.  Effect of Fixing Earthquake Depth in ShakeAlert Algorithms on Performance for Intraslab Earthquakes Effect of Fixing Earthquake Depth in ShakeAlert Algorithms on Performance for Intraslab Earthquakes.  Seismological Research Letters931277–287.  10.1785/0220210056
  * Thompson . (2020) thompsonBasinAmplificationEffects2020Thompson, M., Wirth, EA., Frankel, AD., Renate Hartog, J. Vidale, JE.  2020\.  Basin Amplification Effects in the Puget Lowland, Washington, from Strong-motion Recordings and 3D Simulations Basin amplification effects in the Puget Lowland, Washington, from strong-motion recordings and 3D simulations.  Bulletin of the Seismological Society of America1102534–555.  10.1785/0120190211
  * Tromp . (2010) trompRealtimeSimulationsGlobal2010Tromp, J., Komatitsch, D., Hjörleifsdóttir, V., Liu, Q., Zhu, H., Peter, D.Hutko, A.  2010\.  Near Real-Time Simulations of Global CMT Earthquakes Near real-time simulations of global CMT earthquakes.  Geophysical Journal International1831381–389.  10.1111/j.1365-246X.2010.04734.x
  * Trugman Shearer (2018) trugmanStrongCorrelationStress2018Trugman, DT. Shearer, PM.  2018\.  Strong Correlation between Stress Drop and Peak Ground Acceleration for Recent M 1–4 Earthquakes in the San Francisco Bay Area Strong correlation between stress drop and peak ground acceleration for recent M 1–4 earthquakes in the San Francisco Bay Area.  Bulletin of the Seismological Society of America1082929–945.  10.1785/0120170245
  * Uphoff . (20172) Uphoff2017Uphoff, C., Rettenberger, S., Bader, M., Madden, E., Ulrich, T., Wollherr, S. Gabriel, AA.  20172\.  Extreme scale multi-physics simulations of the tsunamigenic 2004 Sumatra megathrust earthquake Extreme scale multi-physics simulations of the tsunamigenic 2004 Sumatra megathrust earthquake.  Proceedings of the International Conference for High Performance Computing, networking, Storage and Analysis, SC 2017. Proceedings of the international conference for high performance computing, networking, storage and analysis, sc 2017.  https://doi.org/10.1145/3126908.3126948
  * Uphoff . (20171) uphoff2017extremeUphoff, C., Rettenberger, S., Bader, M., Madden, EH., Ulrich, T., Wollherr, S. Gabriel, AA.  20171\.  Extreme scale multi-physics simulations of the tsunamigenic 2004 sumatra megathrust earthquake Extreme scale multi-physics simulations of the tsunamigenic 2004 sumatra megathrust earthquake.  Proceedings of the international conference for high performance computing, networking, storage and analysis Proceedings of the international conference for high performance computing, networking, storage and analysis ( 1–16).
  * van Driel . (2015) vandrielInstaseisInstantGlobal2015van Driel, M., Krischer, L., Stähler, SC., Hosseini, K. Nissen-Meyer, T.  2015\.  Instaseis: Instant Global Seismograms Based on a Broadband Waveform Database Instaseis: Instant global seismograms based on a broadband waveform database.  Solid Earth62701–717.
  * Virtanen . (2020) virtanenSciPyFundamentalAlgorithms2020Virtanen, P., Gommers, R., Oliphant, TE., Haberland, M., Reddy, T., Cournapeau, D.Bright, J.  2020\.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python SciPy 1.0: Fundamental algorithms for scientific computing in Python.  Nature methods173261–272.
  * Wald . (2008) waldQuantifyingQualifyingUSGS2008Wald, DJ., Lin, KW. Quitoriano, V.  2008\.  Quantifying and Qualifying USGS ShakeMap Uncertainty Quantifying and qualifying USGS ShakeMap uncertainty.  US Geological Survey Reston, VA.
  * Wald . (2011) waldUSGSDidYou2011Wald, DJ., Quitoriano, V., Worden, CB., Hopper, M. Dewey, JW.  2011\.  USGS “Did You Feel It?” Internet-based Macroseismic Intensity Maps USGS “Did You Feel It?” Internet-based macroseismic intensity maps.  Annals of Geophysics546.  10.4401/ag-5354
  * Wang Zhan (2020) wangMoving1D3D2020Wang, X. Zhan, Z.  2020\.  Moving from 1-D to 3-D Velocity Model: Automated Waveform-Based Earthquake Moment Tensor Inversion in the Los Angeles Region Moving from 1-D to 3-D velocity model: Automated waveform-based earthquake moment tensor inversion in the Los Angeles region.  Geophysical Journal International2201218–234.
  * Willcox Peraire (2002) willcoxBalancedModelReduction2002Willcox, K. Peraire, J.  2002\.  Balanced Model Reduction via the Proper Orthogonal Decomposition Balanced model reduction via the proper orthogonal decomposition.  AIAA journal40112323–2330.
  * Withers . (2020) withersMachineLearningApproach2020aWithers, KB., Moschetti, MP. Thompson, EM.  2020\.  A Machine Learning Approach to Developing Ground Motion Models from Simulated Ground Motions A machine learning approach to developing ground motion models from simulated ground motions.  Geophysical Research Letters476e2019GL086690.
  * Withers . (2018) withersValidationDeterministicBroadband2018Withers, KB., Olsen, KB., Shi, Z. Day, SM.  2018\.  Validation of Deterministic Broadband Ground Motion and Variability from Dynamic Rupture Simulations of Buried Thrust Earthquakes Validation of deterministic broadband ground motion and variability from dynamic rupture simulations of buried thrust earthquakes.  Bulletin of the Seismological Society of America1091212–228.  10.1785/0120180005
  * Worden . (2018) wordenSpatialSpectralInterpolation2018aWorden, CB., Thompson, EM., Baker, JW., Bradley, BA., Luco, N. Wald, DJ.  2018\.  Spatial and Spectral Interpolation of Ground-motion Intensity Measure Observations Spatial and spectral interpolation of ground-motion intensity measure observations.  Bulletin of the Seismological Society of America1082866–875.  10.1785/0120170201
  * Worden Wald (2016) wordenShakeMapManualShakeMap2016Worden, CB. Wald, DJ.  2016\.  ShakeMap 4 Manual — ShakeMap Documentation Documentation. ShakeMap 4 Manual — ShakeMap Documentation documentation.  http://usgs.github.io/shakemap/.
  * Xiao . (2015) xiaoNonintrusiveReducedorderModelling2015Xiao, D., Fang, F., Pain, C. Hu, G.  2015\.  Non-Intrusive Reduced-Order Modelling of the Navier–Stokes Equations Based on RBF Interpolation Non-intrusive reduced-order modelling of the Navier–Stokes equations based on RBF interpolation.  International Journal for Numerical Methods in Fluids7911580–595.
  * Yang . (2023) yang2023rapidYang, Y., Gao, AF., Azizzadenesheli, K., Clayton, RW. Ross, ZE.  2023\.  Rapid Seismic Waveform Modeling and Inversion with Neural Operators Rapid seismic waveform modeling and inversion with neural operators.  IEEE Transactions on Geoscience and Remote Sensing.
  * Yang . (2021) yangSeismicWavePropagation2021aYang, Y., Gao, AF., Castellanos, JC., Ross, ZE., Azizzadenesheli, K. Clayton, RW.  2021\.  Seismic Wave Propagation and Inversion with Neural Operators Seismic wave propagation and inversion with neural operators.  The Seismic Record13126–134.