See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/304381911

Demonstration of the Cascadia G-FAST Geodetic Earthquake Early Warning

System for the Nisqually, Washington, Earthquake

Article in Seismological Research Letters · June 2016

DOI: 10.1785/0220150255

CITATIONS READS

37 210

11 authors, including:

Brendan W. Crowell Paul Bodin

University of Washington Seattle University of Washington Seattle
59 PUBLICATIONS 1,008 CITATIONS 123 PUBLICATIONS 3,884 CITATIONS

SEE PROFILE SEE PROFILE

John E Vidale Victor Santillan

University of Washington Seattle Central Washington University
228 PUBLICATIONS 9,207 CITATIONS 29 PUBLICATIONS 348 CITATIONS

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Co-authored paper View project

From Boutique to Wholesale Seismic Monitoring: Performance Evaluation Tools to Prepare a Traditional Regional Seismic Network for Earthquake Early Warning View
project

All content following this page was uploaded by Brendan W. Crowell on 10 April 2018.

The user has requested enhancement of the downloaded file.

Demonstration of the Cascadia G-FAST
Geodetic Earthquake Early Warning System
for the Nisqually, Washington, Earthquake
by Brendan W. Crowell, David A. Schmidt, Paul Bodin, John E. Vidale,
Joan Gomberg, J. Renate Hartog, Victor C. Kress, Timothy I. Melbourne,
Marcelo Santillan, Sarah E. Minson, and Dylan G. Jamison
ABSTRACT
A prototype earthquake early warning (EEW) system is cur- imize long-period drifts due to sensor rotations and tilts during
rently in development in the Pacific Northwest. We have taken integration to velocity and displacement (Melgar, Bock, et al.,
a two-stage approach to EEW: (1) detection and initial charac- 2013). This filtering reduces the low-frequency content of the
terization using strong-motion data with the Earthquake Alarm earthquake recording, which is critical to characterize large
Systems (ElarmS) seismic early warning package and (2) the trig- earthquakes.
gering of geodetic modeling modules using Global Navigation Geodetic observations provide an additional constraint that
Satellite Systems data that help provide robust estimates of does not saturate for large-magnitude events. Global Positioning
large-magnitude earthquakes. In this article we demonstrate System (GPS) displacement waveforms on their own have been
the performance of the latter, the Geodetic First Approximation shown to be useful for near-real-time magnitude determination
of Size and Time (G-FAST) geodetic early warning system, using using peak ground displacement (PGD) scaling relationships
simulated displacements for the 2001 M w 6.8 Nisqually earth- (Crowell et al., 2013; Melgar et al., 2015), rapid Centroid
quake. We test the timing and performance of the two G-FAST Moment Tensor (CMT) determination for size and orientation
source characterization modules, peak ground displacement scal- (Melgar et al., 2012; Melgar, Crowell, et al., 2013; O'Toole et al.,
ing, and Centroid Moment Tensor-driven finite-fault-slip mod- 2013), and finite-fault methods that use rapidly computed co-
eling under ideal, latent, noisy, and incomplete data conditions. seismic offsets to compute slip on the fault (Crowell et al., 2009,
We show good agreement between source parameters computed 2012; Allen and Ziv, 2011; Ohta et al., 2012; Wright et al.,
by G-FAST with previously published and postprocessed seismic 2012; Böse et al., 2013; Colombelli et al., 2013; Minson et al.,
and geodetic results for all test cases and modeling modules, and 2014; Grapenthin et al., 2014a). PGD scaling has been shown to
we discuss the challenges with integration into the U.S. Geologi- be slower than P-wave-based methods because (1) large earth-
cal Survey’s ShakeAlert EEW system. quakes take time to reach their peak amplitudes and (2) the
S waves that determine PGD travel little more than half the
speed of P waves. In the case of a large offshore event like the
INTRODUCTION 2011 M w 9.0 Tohoku-Oki earthquake, an initial PGD-based
magnitude estimate of M 8.5 would be available 50 s following
Earthquake early warning (EEW) systems provide seconds to mi- nucleation, which still would have been useful for warning large
nutes of advanced warning after rupture initiation and prior to population centers like Tokyo (Melgar et al., 2015). Methods
the arrival of strong ground shaking at a location. The method- that rely on coseismic offsets can have an impact on traditional
ology currently employed uses seismic networks to quickly char- EEW only in extreme cases (i.e., large earthquakes far from pop-
acterize the magnitude and epicenter from P waves at the ulation centers with S-wave travel times on the order of a few
stations nearest to the epicenter and sends an alert over a broad minutes or preferential geometrical alignment of stations), but
region before the S wave arrives. A common problem with seis- their major utility is in proper characterization of earthquake
mically derived magnitudes for EEW is the saturation of the sig- impact for tsunami early warning and hazard response (Blewitt
nal at large magnitudes (M > 7; Brown et al., 2011). One cause et al., 2006; Ohta et al., 2012). In addition, Crowell et al. (2013)
of saturation is the limited P-wave time window available for have shown that when strong-motion accelerations and GPS dis-
analysis before the arrival of the S wave at nearby stations, which placements are combined using a multirate Kalman filter (Bock
can be shorter than the duration of the rupture source for large et al., 2011), the impact of magnitude saturation is greatly di-
events. Another culprit in magnitude saturation is the high-pass minished, even using only the first 5 s of data after the P-wave
filtering of strong-motion records, which is performed to min- arrival. This method requires collocated GPS and strong-

930 Seismological Research Letters Volume 87, Number 4 July/August 2016 doi: 10.1785/0220150255
motion instruments, which is currently rare in western North currently under development by the U.S. Geological Survey
America. (USGS) with university partners (Given et al., 2014).
The subduction environment of Cascadia motivates the
need for a joint seismic and geodetic EEW system to rapidly
characterize all possible earthquakes within the region. Poten- JOINT EARTHQUAKE EARLY WARNING SYSTEM
tial events include large megathrust events and outer-rise
events anywhere from the Mendocino triple junction to Van- The current ShakeAlert system uses three different seismic al-
couver Island, shallow crustal earthquakes in the Seattle and gorithms for estimating location, magnitude, and origin time
Portland metropolitan areas, and deep events within the slab. (OT): ElarmS (Kuyuk et al., 2014), Virtual Seismologist (VS,
At the University of Washington (UW), development of the Cua and Heaton, 2007), and OnSite (Böse et al., 2009). The
Geodetic First Approximation of Size and Time (G-FAST) California system runs all three algorithms, whereas the Pacific
system has been a high priority for operational EEW in the Northwest system currently only uses ElarmS. Several GPS-
Cascadia region. based algorithms that utilize rapidly computed coseismic offsets
G-FAST continuously receives real-time processed GPS time are currently under development, including G-larmS (Grape-
series from the Pacific Northwest Geodetic Array (PANGA) and nthin et al., 2014a) and BEFORES (Minson et al., 2014). In
maintains a local data buffer. G-FAST receives event triggers this article, we document an additional geodetic module under
from a seismic detection module, such as Earthquake Alarm development, G-FAST. What distinguishes the G-FAST ap-
Systems (ElarmS), and starts with an estimate of the location, proach is its combination of different types of analyses. It first
timing, and size of the earthquake. G-FAST first estimates mag- estimates source depth and magnitude from PGD. This is sim-
nitude and depth from PGD scaling. It then invokes a model- ilar to GPSlip (Böse et al., 2013), which also determines source
ing suite to estimate a CMT and finite-fault parameters. strength from PGD, but with the addition of solving for source
To explore the performance of the system, we imple- depth as well, giving G-FAST an advantage in subduction zone
mented a test system that reads in synthetic data and xml mes- environments where earthquakes occur over a wide range of
sages from a seismic detection module and outputs the depths. G-FAST then uses a geodetically derived focal mecha-
nism to determine the fault orientation, builds a discretized
information in a simulated real-time mode, leaving the
fault plane with that orientation, and then inverts for the static
back-end modeling modules untouched. The simulated system
slip distribution on that fault plane. This final slip model is as
allows us to vary the latency, data completeness, and noise to
complex as the BEFORES geodetic EEW algorithm (Minson
test the robustness of magnitude, timing, and slip estimates
et al., 2014), which rapidly updates both the fault orientation
from G-FAST.
and the spatial distribution of accumulated slip on that fault
In this article, we first give an overview of G-FAST and the
plane as the rupture evolves, and is more complex than the G-
synthetic test system. Then, we demonstrate the performance
larmS geodetic EEW source model (Grapenthin et al., 2014a),
during a simulation of the 28 February 2001 M w 6.8 Nisqually
which only allows for along-strike variations in slip amplitude.
earthquake. The Nisqually earthquake was a deep intraslab
Thus, the G-FAST slip model is probably better suited than G-
event located at the southern end of Puget Sound, roughly larmS to subduction zone earthquakes that can have significant
50 km deep, and caused damage costing several billion dollars along-dip variation in slip. The G-FAST approach is also sub-
in the Seattle area (Ichinose et al., 2004). The Nisqually earth- stantially faster than G-larmS because we can obtain the PGD-
quake is a good test case for G-FAST, because intraplate events determined magnitude and source depth before the wavefield
have a higher probability of occurrence in the Pacific North- has converged to the static state, and it is probably comparable
west than a megathrust event (∼30–50 year recurrence, with in speed to BEFORES, although direct head-to-head perfor-
the prior comparable 1949 Olympia and 1965 Seattle–Tacoma mance testing between the various geodetic EEW algorithms
events, Ichinose et al., 2004), and it caused shaking over a wide has yet to be done. The flowchart of operations of G-FAST is
region (modified Mercalli intensity [MMI] VI–VII through- shown in Figure 1. The system consists of two independent
out the Puget Lowlands, see Data and Resources), and was re- packages: (1) a data aggregator and buffer that runs continu-
corded with fairly low signal-to-noise surface displacements. ously and (2) a triggered modeling suite that activates upon
The components of G-FAST have been independently tested receiving an alert from the seismic warning system (currently
and validated on much larger earthquakes elsewhere in Japan, ElarmS). High-rate real-time GPS positions are continuously
Chile, Indonesia, and southern California (Crowell et al., received at the Pacific Northwest Seismic Network (PNSN)
2012, 2013; Melgar et al., 2012, 2015); so the motivation for from PANGA through the JSON protocol (see Data and
testing G-FAST on the Nisqually earthquake is to determine Resources). The PANGA solutions are integer ambiguity-
the performance and resolution toward the lower end of de- resolved precise point positions (Zumberge et al., 1997) esti-
tectability and to ascertain statistics on the range of possible mated at 1 s epochs within the ITRF2008 reference frame (Al-
solutions by varying the latency, noise, and data completeness. tamimi et al., 2012) using satellite orbit and clock corrections
Finally, we discuss how to best improve the EEW system based provided by the International Global Navigation Satellite Sys-
on the simulation results and challenges associated with inte- tems Service. Station positions are estimated independently
gration into the ShakeAlert system, which is the EEW system and do not depend on a fixed reference station or network.

Seismological Research Letters Volume 87, Number 4 July/August 2016 931

▴ Figure 1. Flowchart of proposed joint early warning system with Geodetic First Approximation of Size and Time (G-FAST) and seismic
algorithms (Earthquake Alarm Systems [ElarmS]). The color version of this figure is available only in the electronic edition.

We maintain a continuous 5-min data buffer as well as archive epicentral region (Bustin et al., 2004), although none of these
data at UW. were capable of recording high-rate data (1 Hz or greater). Be-
ShakeAlert’s ElarmS software continuously scans strong cause of this, we supplement the strong-motion recordings with
motion and broadband data from the PNSN and other nearby a set of synthetic displacement waveforms at the locations of
seismic networks for possible P-wave arrivals. When ElarmS the 26 strong-motion stations, computed up to 5 Hz using the
identifies a possible event at four seismic stations, it determines frequency–wavenumber (f
k) integration method (Zhu and
an epicenter, OT, and magnitude with the first 4 s or less of Rivera, 2002). The method sums the waves radiated from an
data after the P-wave arrival at each station. ElarmS computes array of point sources on a predefined fault to model wave-
peak displacements 10 times per second and updates magni- forms. The fault we use is centered on the PNSN hypocenter
tude and epicenter as more stations report arrivals. When an location (47.149° N, 122.727° W, 51.9 km deep) with a strike
M 3.0 or greater event is reported, the geodetic modeling mod- of 350°, dip of 70°, rake of −90°, length of 23 km, width of
ules are initiated using the OT and epicenter from ElarmS. This 10 km, and magnitude 6.8 (Bustin et al., 2004). For the seismic-
magnitude is well below the noise threshold of detection with velocity structure, we use the standard layered half-space model
real-time GPS (M ∼ 6), but we choose to trigger on all events that is used by the PNSN to locate earthquakes in Puget Sound.
reported by the seismic warning system. The PGD module es- The shear modulus is assumed to be 66 MPa in the seismic
timates magnitude and depth once four GPS stations are within source region (Bustin et al., 2004).
a 3
km=s travel-time mask and updates these estimates every The amplitudes and timing of the synthetic velocity wave-
second. The CMT-derived finite-fault inversion waits until forms are in qualitative general agreement with the recorded
static offsets (1
km=s travel-time mask with 10 s of averaging) seismic waveforms, although precise matches were not achieved
are in hand to compute an estimate of fault orientation, mag- nor needed for our purposes. The peak velocity root mean
nitude, depth, and slip along the fault surface. The source mod- square (rms) difference between the synthetics and the data was
els from both modules will be fed into a joint seismic/geodetic 0:05 m=s across 69 channels, and 60% of the synthetic channels
ShakeAlert DecisionModule that, although yet to be devel- achieved peak velocity within 5 s of the data. For our synthetic
oped, will combine the information from all of the contribut- source model, a constant fault slip of 1.4 m is assigned across
ing EEW algorithms to produce a unified shaking forecast. This the fault plane with a rupture front that propagates from the
shaking forecast would then be sent to end users. hypocenter at 3:2 km=s. We did not attempt to account for any
additional complexity of the rupture, which might better em-
DATA FROM THE 2001 NISQUALLY EARTHQUAKE ulate the observed waveforms. For G-FAST, a simpler source is
sufficient because the largest impact on source estimation will
At the time of the Nisqually earthquake, the PNSN consisted of be the amplitudes of the displacement waveforms, the timing of
26 strong-motion stations in the region (Fig. 2). The PANGA the peak displacement, and the time in which the static offset is
network consisted of ∼15 permanent GPS stations within the fully emplaced. The predicted synthetic displacements agree

932 Seismological Research Letters Volume 87, Number 4 July/August 2016

 48˚N

UW.LAWT SEAT

47˚N UW.RWW/SATS

km

0 50

1 cm
46˚N gray − FK coseismic

black − Bustin et al. [2004] coseismic

124˚W 123˚W 122˚W 121˚W

▴ Figure 2. The gray vectors indicate the location of the strong-motion stations used in this study as well as their associated
coseismic displacements. The black vectors are coseismic displacements computed from Bustin et al. (2004). The star indicates the
location of the epicenter. Note that three of the strong-motion stations are outside of the map. The color version of this figure is available
only in the electronic edition.

well with the observed static displacements from GPS (Bustin GEODETIC MODELING MODULES
et al., 2004). Figure 2 shows the coseismic displacements from
the synthetics (gray arrows) and the observed offsets deter- Peak Ground Displacement
mined by Bustin et al. (2004) (black arrows). Although there Crowell et al. (2013) presented the scaling of PGD as a func-
were not many GPS stations at the time, GPS station SATS is tion of hypocentral distance for earthquakes between M 5.4
collocated with strong-motion station UW.RWW, and GPS sta- and 9.0. They found the following relationship:
tion SEAT is 6.1 km from strong-motion station UW.LAWT.
log 10 PGD  A  BM w  CM w log 10 r hyp ; 1
The three-component rms differences between simulated and
EQ-TARGET;temp:intralink-;df1;323;166

measured coseismic displacements at these two stations are 2.82 in which r hyp is the hypocentral distance in kilometers, M w is
and 1.35 mm, respectively, and other GPS stations in the region the moment magnitude, PGD is the maximum of the Euclidian
show a similar pattern of deformation to our synthetic coseismic norm of the three components of displacement (north, east,
displacements. We only use simulated displacements from the and vertical) in centimeters, and A, B, and C are the constants
subset of 26 operational strong-motion stations in 2001 because solved through a regression of previous data. The following
ground motions were validated only at those stations. inverse problem is set up to solve for magnitude:

Seismological Research Letters Volume 87, Number 4 July/August 2016 933

 (a) 10 2
(b) 100 (c)
7.3
Correct Depth

80 7.2

Variance Reduction (%)

Shallower Depths Deeper Depths

Magnitude
PGD (cm)

7.1

101 60

7.0

40
6.9

101 102 0 20 40 60 80 100 0 20 40 60 80 100

Hypocentral Distance (km) Depth (km) Depth (km)

▴ Figure 3. (a) The predicted peak ground displacement (PGD) for an M 7.0 earthquake at a depth of 25 km (correct depth in figure and
black circles). The gray squares and white triangles show how the curve would change for different assumed depths (5 km for shallower,
45 km for deeper). (b) The variance reduction (VR) as a function of depth for the PGD regression for the same earthquake. In this example,
at 25 km the VR is 100%, meaning that the model (log–log linearity) fits the data perfectly. (c) The inverted magnitude as a function of depth.

G M w  b;
EQ-TARGET;temp:intralink-;df2;40;505 2 et al. (2015) using an expanded data set of 10 earthquakes re-
corded by GPS alone. The recalibration of the Crowell et al.
2 3 2 3 (2013) regression uses only seismogeodetic (collocated GPS and
B  C log 10 r hyp;1  log 10 PGD1  − A strong-motion stations) data and hence represents a minimal
6 .. 7 6 .. 7
G4
EQ-TARGET;temp:intralink-;df3;40;488

. 5; b4 . 5: noise solution. However, in general, all the regression solutions

B  C log 10 r hyp;n  log 10 PGDn  − A predict similar magnitudes within the bounds of uncertainty of
the method (0:3 magnitude units).
3 Finally, we introduce a grid search for the earthquake
For G-FAST, we make a few key operational changes from depth. ElarmS assumes a depth of 8 km for all events, because
Crowell et al. (2013). First, we introduce a travel-time mask that is the average seismogenic depth in California. But the
PGD algorithm requires proper depth characterization, given
of 3 km=s given the earthquake OTs from ElarmS that ignores
all stations outside of the travel-time mask at a given time. Sec- that the depth affects the PGD pattern. It is also important
ond, no magnitudes are estimated for fewer than four stations. to discriminate between shallow crustal earthquakes and deeper
Third, we utilize exponential distance weighting of the form: events that are located along the plate interface or within the
slab. The grid-search method simply computes the magnitude
  at 1 km intervals between 0 and 100 km depth and chooses the
r 2epi;i
wi  exp − 2 ; 4 depth that maximizes the variance reduction (VR) defined by
8r epi;min
EQ-TARGET;temp:intralink-;df4;40;335

 
in which r epi;i is the epicentral distance in kilometers of the ith jjb − GM w jj
VR  1 − × 100: 6
station and r epi;min is the epicentral distance of the closest sta- jjbjj
EQ-TARGET;temp:intralink-;df6;311;301

tion. The distance weighting in equation (4) is a function of

the epicentral distance, so as not to bias the depth grid search We demonstrate this depth-dependent grid-search ap-
discussed later toward shallower solutions. The factor of 8 in proach for a hypothetical M 7.0 earthquake at 25 km depth
the denominator is arbitrarily chosen through trial and error to (Fig. 3). We use the scaling relationship of equation (1) to find
give relatively high weight to many close stations before drop- the PGD at a number of hypocentral distances, given a 25-km
ping off exponentially. The magnitude is then found by solving depth. Then, we recompute the hypocentral distances based on
different depths (between 1 and 100 km), invert for the mag-
W GM w  Wb;
EQ-TARGET;temp:intralink-;df5;40;191 5 nitude, and compute the VR of the inversion. In a perfect case,
the PGD versus distance plot will be linear in log–log space,
in which W is a diagonal matrix of station weights wi . We recali- shown as correct depth in Figure 3a. For shallower and deeper
brate the regression constants of Crowell et al. (2013) with the depths, the PGD–distance curve (Fig. 3a) is no longer linear on
inclusion of the distance weight matrix using data for the a log–log plot, so a straight line will never fit the data perfectly;
Tohoku-Oki, Tokachi-Oki, and El Mayor–Cucapah earth- this is seen in Figure 3b where the 25 km depth yields a perfect
quakes. We find A  −6:687, B  1:500, and C  −0:214, 100% VR and all other depths fit the model worse. Figure 3c
with an improved magnitude uncertainty of 0.17 magnitude shows how the assumption of depth will change the magnitude
units. These coefficients are different from an analysis by Melgar estimate. For this hypothetical earthquake, changing the source

934 Seismological Research Letters Volume 87, Number 4 July/August 2016

depth for 1–100 km will change the magnitude by 0.4 mag- with the Green’s functions, G ipq;n , being defined by equation (7).
nitude units. Although this magnitude change is not large, the The displacements ui;n are the static offsets on each of the three
prediction of the strength of expected ground shaking is sen- directional components, which are computed by taking the aver-
sitive to different source depths. age of the first 10 s of data that arrive after a travel-time mask of
1 km=s is used. Note, displacements are with respect to the sta-
tion position at the OT of the earthquake. This travel-time mask
CMT-Driven Slip Modeling on a Finite Fault is overly conservative, but we want to minimize inclusion of any
Computing the CMT is important, because it allows us to de- dynamic motions that may contaminate the static offset measure-
termine the fault orientation and location as well as magnitude, ments and do not want to rely on some other metric to deter-
which have important implications for expected ground mo- mine if the solution is stable. The displacements and Green’s
tions and tsunami potential. Melgar et al. (2012) showed how functions are both rotated into the radial, transverse, and vertical
to compute the CMT using rapidly computed static offsets in a directions, and the moment tensor is decomposed into the main
layered half-space. For efficiency and simplicity, we chose to and auxiliary fault planes using the Python ObsPy package (see
solve the problem in a homogeneous half-space using the static Data and Resources). Rather than performing a grid search for
displacement field analytical solution for the moment tensor the centroid location, we use the epicentral location from ElarmS
from Hashima et al. (2008) and the references therein. The gen- and only perform a grid search for the depth using the same VR
eral solution for the displacement field ui x for a given moment maximization scheme as was done for PGD.
tensor M pq at a point x from a source at point ξ is given by A slip inversion on a finite fault provides a more accurate
1 characterization of the event when compared to a point source,
ui x  γ3ζi ζp ζq − ζ i δpq − ζ p δqi − ζ q δip  especially given a potential large megathrust event in Cascadia.
8πμR2
EQ-TARGET;temp:intralink-;df7;52;541

Melgar, Crowell, et al. (2013) outlined the difficulties of using

 2ζq δip M pq ; 7 a simple point-source approximation for computing the CMT
from static offsets for the Tohoku-Oki earthquake. In that case,
in which i, p, and q are the three Cartesian directions, they found using a fixed hypocenter point source would lead to
γ  3K  μ=3K  4μ, μ and K are the rigidity and bulk an accurate moment tensor solution; however, the model fit
modulus, respectively, δip is the Kronecker delta, R is the source– would be poor (VR < 50%) and could not be trusted in real
receiver distance defined by time. A slip inversion on a finite fault does not have these issues
q if the fault plane is large enough, is in the correct region, and
R
EQ-TARGET;temp:intralink-;df8;52;434 x1 − ξ1 2  x2 − ξ2 2  x3 − ξ3 2 ; 8 has reasonable strike and dip angles. Grapenthin et al. (2014a)
investigated the slip inversion sensitivity to location and fault
orientation errors. Although a general rule-of-thumb is diffi-
and the scaled source–receiver distance is
cult here, because this is highly dependent on the earthquake
x i − ξi source and network geometry, they found that variations in lo-
ζi  : 9 cation, dip, and strike become less important for deeper events,
R
EQ-TARGET;temp:intralink-;df9;52;379

and for shallower events the orientations should be within 5°.

The following inverse problem for n stations is set up and solved Misplacement of the center of the fault can impact the magni-
with linear least squares tude error by up to 0.5 magnitude units within 20 km.
For G-FAST’s slip inversion, we use the method of Crowell
2 3 2 3 et al. (2012) where the fault geometry is defined by the CMT,
u1;1 G 111;1 G 121;1 G 131;1 G 221;1 G 231;1 G 331;1
EQ-TARGET;temp:intralink-;df10;52;311

6u 7 6G and the Green’s functions are prescribed by Okada’s formu-

6 2;1 7 6 112;1 G 122;1 G 132;1 G 222;1 G 232;1 G 332;1 7
7 lation (Okada, 1985) in a homogeneous half-space. The center
6 7 6 7
6 u3;1 7 6 G 113;1 G 123;1 G 133;1 G 223;1 G 233;1 G 333;1 7 of the fault plane is defined by the ElarmS epicenter and the
6 7 6 7
6 . 7 6 . .. 7 depth computed from the CMT. The along-strike and along-
6 . 76 . .. .. .. .. 7
6 . 7 6 . . . . . . 7 dip dimensions of the fault are defined by scaling relationships
6 7 6 7
6u 7 6G G 121;n G 131;n G 221;n G 231;n G 331;n 7 from Dreger and Kaverina (2000), based on the CMT magni-
6 1;n 7 6 111;n 7
6 7 6 7 tude. We are not concerned that the CMT magnitude may be
4 u2;n 5 4 G 112;n G 122;n G 132;n G 222;n G 232;n G 332;n 5
an overestimate, because this will only make the potential slip
u3;n G 113;n G 123;n G 133;n G 223;n G 233;n G 333;n surface larger; the inversion does not prescribe slip in areas if the
2 3
M 11 data misfit does not call for it. The CMT-computed depth also
6 M 12 7 will not greatly impact the result, because it defines the center of
6 7 the fault and the along-dip dimension will cover the majority of
6 7
6 M 13 7 the possible seismogenic zone. Of the two fault planes from the
×6
6 M 7;
7 10
6 22 7 CMT, we pick the one that minimizes the misfit as the final
6 7 solution. Laplacian regularization with a generalized smoothing
4 M 23 5
parameter described in Crowell et al. (2012) is used. The static
M 33 offsets used are the same as for the CMT inversion.

Seismological Research Letters Volume 87, Number 4 July/August 2016 935

SIMULATIONS RESULTS
We first run G-FAST assuming no data latency and no data Peak Ground Displacement
noise to find the ideal source characterization and timing with The evolution of PGD magnitude and depth estimates over
the synthetic displacements. The ElarmS location, timing, and time is shown in Figure 4 for the ideal case as well as the four
magnitude are found using the actual strong-motion data be- simulations for the Nisqually earthquake. In the ideal case, the
cause the short-term average/long-term average algorithm for first magnitude estimates from PGD are available 17 s after the
detecting P-wave arrivals requires some level of noise to operate OT, trailing ElarmS by about 4 s. The initial PGD magnitude
efficiently. estimate is M 6.47 (at 17 s) and M 6.43 for ElarmS (at 12.8 s).
To obtain realistic operational system performance, we The arrival of strong shaking in Seattle, 60 km northwest of the
run four simulations and for each conduct 1000 trials of G- epicenter, is at 23 s after OT, so ElarmS would provide a 10 s
FAST. The four simulations test the impacts of (1) randomly warning to Seattle and an updated warning from PGD is avail-
generated latencies, (2) high-rate GPS station noise, (3) data able 6 s prior to strong shaking. For both ElarmS and PGD, the
dropouts, and (4) latency + noise + dropouts. We perform magnitude estimates start out small and increase quickly over
these simulations to test the stability and robustness of the the first ∼10 s to their final stable magnitude estimate. In the
source solutions, obtain realistic timing and to see how each ideal case, no PGD estimate is provided to locations within
factor impacts the results for future improvements. We gener- 51 km of the epicenter. The magnitude estimates over time
ate integer latencies from a Poisson distribution with a mean of do not vary much for the ideal case; the full range of magnitude
6 s. Data latency from the PANGA/Plate Boundary Observa- estimates is 0.2 magnitude units. The stable magnitude esti-
tory (PBO) network is generally much better than this and is mate (> 30 s) from PGD (M 6:7  0:3) and from ElarmS
improving (phase and range distribution < 1 s, processing and (M 6.9) are both close to the true magnitude of 6.8.
distribution < 3 s); however, we choose a conservative estimate. The depth estimate in the ideal case starts out shallow but
We generate the station noise using the power spectrum quickly converges to the final solution of 52 km, exactly the
for 40-km relative GPS positions from Genrich and Bock input model depth (51.9 km) and close to the Global CMT
(2006), a combination of white (f 0 ), flicker (1=f ), and ran- solution of 46.8 km. ElarmS assumes a depth of 8 km because
dom-walk (1=f 2 ) noise. We choose to simulate noise instead it was designed to represent the average seismogenic depth in
of superimposing recorded noise time series from PANGA to California; this has an impact on the prediction of travel times
have control over the range of potential noise sources. We of strong shaking as well as estimating ground motions in
compared the power spectra of the 10 most complete time ShakeAlert. Assuming a 1=r attenuation, the 8 km depth at an
series recorded at the PNSN in real time from PANGA on epicentral distance of 20 km would overestimate strong shak-
15 September 2015 with the power spectra of Genrich and ing by a factor of 2.5 versus the 50 km depth given the same
Bock (2006). At periods less than 10 s, the PANGA real-time magnitude earthquake. The assumption of a shallower depth
solutions matched perfectly with the Genrich and Bock (2006) would lead to a smaller PGD magnitude estimate, however, as
power spectra. At periods between 5 min and 10 s, the real-time demonstrated by Figure 3c. We plan on modifying ElarmS in
power spectra are systematically lower than the simulation power the Pacific Northwest to accommodate a depth grid search,
spectra we use, indicating that the simulations in this article re- however, dealing with other aspects of tuning ElarmS for the
present a worst-case scenario, especially at long periods. region has taken precedence (Hartog et al., 2016; see Data and
The real-time data completeness rate for PBO in the Resources).
Cascadia region is generally better than 95% (D. Mencin, All four simulations produce stable estimates of magnitude
personal comm., 2015, UNAVCO). Most data dropouts are after about 30 s, although the range of magnitude estimates
due to a few stations with less than ideal telemetry, but for the prior to 30 s is less than a magnitude unit, and the range of
sake of argument, we assume a data return rate of 85% for all solutions quickly narrows to ∼0:3 magnitude units. Evidence
stations during our simulations. For a large earthquake, we are of this tight distribution of magnitude estimates is shown in
uncertain as to the telemetry robustness, and therefore we want Figure 5, which contains histograms of the simulations at
to consider a worst-case scenario. For these simulations, we 30 s after OT. The impact of latency is minimal on the stability
simply remove 15% of the data points at random each trial. of magnitude, and its impact is concentrated toward the begin-
We do not consider spatial or temporal correlation of dropouts ning of the earthquake, although latency is the only parameter
even though a significant percentage of dropouts will be due that impacts the warning time. The average first-alert time in
to the temporary failure of a single telemetry path that may the latency simulations is 21.9 s after OT, which matches the
handle several stations. Even without explicitly considering 5 s Poissonian distributed latency used in the simulations. By
spatially correlated dropouts, removing 15% of the data (when about 70 s after OT, the latency simulations are exactly the
data completeness is generally much better than 95%) from same as the ideal case. This is not surprising, because latency
each station at random will explore the range of possible sol- will only alter the order in which data come in and are utilized;
utions, and the source parameters should be robustly estimated after some time, all the important data (i.e., the recording of
with 1000 trials. peak ground motions) will be available. Latency does, however,

936 Seismological Research Letters Volume 87, Number 4 July/August 2016

 (a) (e) 100
7.0
ElarmS 90

6.8 80

PGD 70
6.6
60

6.4 50
Latency 40
6.2
30

6.0 20

(b) (f) 100

7.0
90

6.8 80

70
6.6
60

6.4 50
Noise 40
6.2
30

Depth (km)
Magnitude

6.0 20

100

(c)7.0 (g) 90

6.8 80

70
6.6
60

6.4 50
Dropouts 40
6.2
30

6.0 20

(d) (h) 100

7.0
90

6.8 80

70
6.6
60

6.4 50
Latency+Noise+Dropouts 40
6.2
30

6.0 20
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
Seconds after OT

▴ Figure 4. (a–d) Magnitude and (e–h) depth as a function of time after origin time (OT) for the PGD simulations. The thick black lines
represent the ideal solution with no noise, latency, or dropouts; the dashed lines on the left plots show the ElarmS magnitude estimates
over time, and the thin gray lines represent different iterations in the four simulations, (a,e) latency, (b,f) noise, (c,g) dropouts, and (d,
h) latency + noise + dropouts.

delay the first-alert time from PGD by about 4–5 s, thus delay- timate for that station by 0.28 magnitude units; however, a
ing the first PGD estimate of magnitude. station with a PGD of 10 cm at 100 km and an additional
The impact of noise is also not surprising. For the most 1 cm of PGD added through noise will only increase the mag-
part, noise causes the estimate of magnitude to be greater than nitude estimate by 0.04 magnitude units. Dropouts have the
in the ideal case, because the time series are just a superposition greatest impact on the magnitude estimates and tend to bias
of the noise and the signal. Many of the stations have dynamic toward lower-magnitude estimates than the ideal case, due to
displacements that are not much greater than the noise level, so the removal of peak ground motions; however, as evidenced by
this low signal-to-noise ratio (SNR) can have a significant im- the histogram in Figure 5c, the vast majority of magnitude es-
pact on the magnitude estimates. For example, a station that timates at 30 s after OT in the simulation fall in a 0.2-magni-
has a PGD of 1 cm at 100 km distance and an additional 1 cm tude band around M 6.6. The simulation looking at latency,
of PGD added through noise will increase the magnitude es- noise, and dropouts has a wider band of possible magnitude

Seismological Research Letters Volume 87, Number 4 July/August 2016 937

(a) 100 (e) 30 s are between 40 and 65 km depth, which encompasses the
previously published depth estimates (Bustin et al., 2004; Ichi-
Latency nose et al., 2004; Kao et al., 2008). Dropouts cause the greatest
impact on the depth estimate, due to the removal of peak dis-
50 placements that would impact the model fit. The distribution
of depth measurements about the ideal solution for dropouts
exhibits a skewing toward shallower depths. This is under-
standable, because the dropouts will impact the PGD measure-
0
6.5 7.0 20 40 60 80 100 ments from the closest stations most, which will reduce the
(b) 100 (f) slope of the PGD–distance curve. However, most of the depth
estimates fall between 30 and 60 km in the case of dropouts.
Noise
The histogram for the case of latency, noise, and dropouts
shows a similar shape to the dropout histogram. This implies
50
that dropouts are the most important factor to consider, fol-
lowed closely by noise. The effect of latency is minimal, and
reducing latency will only have an impact on warning time and
Percent

0 not on the stability of the magnitude and depth estimates.

6.5 7.0 20 40 60 80 100
100
(c) (g)
CMT-Driven Slip Inversion
Dropouts
The CMT simulation results showing the magnitude, depth,
strike, dip, and rake evolution as a function of time are shown
50
in Figure 6. The first CMT (and finite-fault) estimate is avail-
able 38 s after OT for the ideal case, with magnitude M 7.3 and
strikes, dips, and rakes of 190°/319°, 11°/82°, and −40°= − 99°,
0
respectively. After 50 s, the CMT source parameters only vary a
6.5 7.0 20 40 60 80 100 small amount and can be viewed as fully stable after this time.
(d) 100 (h) The ideal strikes, dips, and rakes at 100 s are 216°/344°,
16°/81°, and −38°= − 102°, which is in general agreement with
Latency+Noise+Dropouts
the USGS solution of 172°/347°, 20°/71°, and −85°= − 92°
50
and the input model of 350° strike, 70° dip, and −90° rake.
We find an ideal source depth is 49 km, which is in line with
previous published results placing it at ∼50 km and the input
model depth of 52 km. Our final ideal CMT magnitude is
0 M 7.1, larger than the actual magnitude of 6.8, but understand-
6.5 7.0 20 40 60 80 100
able given the simplifying assumptions of the method (point
Magnitude Depth (km)
source) and the low number and low distribution of stations
being used to compute the CMT.
▴ Figure 5. Histograms of (a–d) magnitude and (e–h) depth for
For the CMT simulations, the impact of latency is similar
the four different PGD simulations at 30 s after OT. (a,e) Latency
to the PGD results, in that the variability in the source param-
simulations, (b,f) noise, (c,g) dropouts, and (d,h) latency + noise +
eters is greater toward the beginning of the recordings, with a
dropouts for magnitude and depth, respectively. The bin width for
similar 5 s delay in the first-alert time. Noise does not lead to a
magnitude is 0.05 magnitude units and for depth is 5 km. The ver-
bias, as was the case for the PGD simulations, but rather leads
tical lines represent the optimal PGD solution at 30 s after OT, mag-
to a fairly even spread about the ideal solution across all source
nitude 6.65, and depth 45 km.
parameters. Dropouts on their own do not impact the CMT
results at all. This is not surprising, because we average 10 s of
estimates of about 0.3 magnitude units (Fig. 4d), but the result data to determine the static offsets, so a few missing data points
is still robust and only slightly larger than the previously pub- will not impact this average when there is no noise. If all data
lished uncertainties of the PGD method (Crowell et al., 2013; points are missing, G-FAST removes the station from further
Melgar et al., 2015). calculations, and any solution variability will be due to station
The depth results have a much greater range of possible distribution changes. Note, this in effect tests the temporally
solutions, but like the magnitude estimates, the histograms correlated dropouts scenario. However, the combined effects of
in Figure 5 show the distributions are much tighter than they noise and dropouts do provide appreciable variability as evi-
appear in Figure 4. Latency, once again, shows that it has the denced by the simulation under all conditions. The CMT sim-
smallest impact on the results and eventually converges to the ulations under all conditions provide insight into the statistical
ideal solution by 70 s. The noise simulations bias the depth uncertainties of the CMT estimation. We compute the stan-
estimates toward deeper depths, although all the simulations at dard deviations of each parameter (magnitude, depth, strike,

938 Seismological Research Letters Volume 87, Number 4 July/August 2016

 Magnitude Depth (km) Strike Dip Rake
(a1) (b1) 100 (c1) 360 (d1) 90 (e1) 0
7.4
90
7.2 300
CMT 80
7.0 240 60 −60
ElarmS 70
6.8 60 180
6.6 50
120 30 −120
6.4 40
Latency 60
6.2 30
6.0 20 0 0 −180
(a2) 7.4 (b2) 100 (c2) 360 (d2) 90 (e2) 0
90
7.2 300
80
7.0 240 60 −60
70
6.8 60 180
6.6 50
120 30 −120
6.4 40
Noise 60
6.2 30
6.0 20 0 0 −180
(a3) (b3) 100 (c3) 360 (d3) 90 (e3) 0
7.4
90
7.2 300
80
7.0 240 60 −60
70
6.8 60 180
6.6 50
120 30 −120
6.4 40
Dropouts 60
6.2 30
6.0 20 0 0 −180
(a4) 7.4 (b4) 100 (c4) 360 (d4) 90 (e4) 0
90
7.2 300
80
7.0 240 60 −60
70
6.8 60 180
6.6 50
120 30 −120
6.4 40
Latency+Noise+Dropouts 60
6.2 30
6.0 20 0 0 −180
30 40 50 60 70 80 90 100 30 40 50 60 70 80 90 100 30 40 50 60 70 80 90 100 30 40 50 60 70 80 90 100 30 40 50 60 70 80 90 100
Seconds after OT

▴ Figure 6. Centroid Moment Tensor (CMT) simulation results for (a1–a4) magnitude, (b1–b4) depth, (c1–c4) strike, (d1–d4) dip, and (e1–
e4) rake. The thick black lines represent the ideal solution with no noise, latency, or dropouts, and the thin gray lines represent different
iterations in the four simulations, latency (row 1), noise (row 2), dropouts (row 3), and latency + noise + dropouts (row 4).

dip, and rake) at 45 s after OT, because only an average of four analysis, the final VRs under the ideal case are 97% and 92%,
stations are used to compute the CMT at that time. We find indicating that both fault models describe the data exception-
errors of 0.1 magnitude units, 17 km depth, 35.5° in strike, ally well, and differences in shaking prediction are minimal.
12.9° in dip, and 22° in rake. By 60 s after OT, these errors We attribute this problem to both the station density and dis-
decrease by at least a factor of 3. tribution, with most stations on the eastern side of the fault, as
The finite-fault magnitude evolution for the four simula- well as to the low SNR of the stations in this study. The slip
tions is shown in Figure 7. Figure 7 shows that the magnitude inversion without noise also selected the auxiliary fault plane as
derived from the finite fault (∼6:7) is much closer to the real the preferred solution.
magnitude, and the different simulations have similar behavior The bias in the finite-fault solution attributed to the
to the CMT simulations for which noise and all three condi- skewed station distribution is analogous to an offshore event
tions cause the greatest variability. The magnitude estimates in which stations are on land toward the east of an event. In
from the finite-fault simulations have less variability than the Figure 8, the current seismic and real-time GPS networks are
CMT simulations. The optimal slip model and GPS fits at 100 s shown by the gray symbols. To test whether or not we would
are shown in Figure 8. At 100 s after OT, only 15 of the sta- expect issues with the fault-plane orientations with our current
tions are within the travel-time mask. The most troubling issue network, we forward modeled coseismic static displacements at
with the slip inversion is the preference for the shallowly dip- the current GPS network stations and the old strong-motion
ping fault plane, which is not thought to be the correct plane network configuration due to slip on the main (steeply dip-
and is not the input plane for the f
k integration. Several ping) fault plane, using an even distribution of 1 m slip (reverse
other studies also struggled with fault-plane ambiguity for this faulting), added random noise, and then inverted the static dis-
event (Ichinose et al., 2004; Kao et al., 2008). Although the placements onto both the main and auxiliary fault planes.
auxiliary fault plane is the preferred plane by the finite-fault Figure 9 shows the difference in the VR of the inversion

Seismological Research Letters Volume 87, Number 4 July/August 2016 939

(a) 7.0 (c) 48˚N

6.9

6.8

6.7

6.6

Latency Dropouts
Magnitude

6.5
40 50 60 70 80 90 100 40 50 60 70 80 90 100
7.0
47˚N
(d)
(b)
6.9

6.8
km
6.7 0 50
1 cm
(black−Data;gray−Model)
6.6 Current PNSN
Noise Latency + Noise + Dropouts Current RTGPS
6.5 46˚N
40 50 60 70 80 90 100 40 50 60 70 80 90 100
Seconds after OT 124˚W 123˚W 122˚W 121˚W

▴ Figure 7. Finite-fault magnitude estimates for the four simula-

45
NE SW
tions; (a) latency, (b) noise, (c) dropouts, and (d) latency + noise +
Depth (km)
dropouts. The thick black lines represent the ideal solution with
no noise, latency, or dropouts, and the thin gray lines represent
different iterations in the four simulations. 50 1m

between the steeply dipping and shallowly dipping fault planes 0 10 20 30 40 50 60

Along Strike (km)
for 1000 iterations. With this even slip distribution, both net-
work configurations choose the steeply dipping plane, but the
▴ Figure 8. The optimal slip model and coseismic fits at 100 s
current (denser) network configuration shows a greater spread
after OT. The optimal slip model is shown in the bottom panel,
in the VRs between the main and auxiliary fault-plane inver-
with the vectors indicating the motion of the hanging wall with
sions than the older network configuration, by an average of
respect to the foot wall. The map view shows the input coseismic
2%. Considering that the VRs are already high, this 2% change,
displacements (black) and modeled displacements (gray). The
although small, is not an insignificant improvement. Although
gray circles show the locations of the current real-time Global
this result is not as satisfying as we would desire, for a medium-
Positioning System (GPS) stations and gray squares are the cur-
sized deep earthquake, the current network configuration still
rent strong-motion stations in the Pacific Northwest Seismic Net-
leads to an improvement in the fault-plane ambiguity, and
work (PNSN). The focal mechanism plot in the map view is the
greater station density would have a similar impact.
ideal CMT computed at 100 s, which the fault plane (black circles)
is centered on. The color version of this figure is available only in
FUTURE DIRECTIONS FOR GEODETIC EEW the electronic edition.

One limitation with G-FAST lies with the lower limit of de- shocks of the 2014 Napa earthquake and showed that G-larmS
tectability. With GPS data alone at regional distances, any will consistently compute M ∼ 6 for small earthquakes due to
event M < 6 will be difficult to measure due to noise. With noise, but the misfit of the models indicate that the magnitude
seismogeodetic (GPS + strong-motion) data, Geng et al. (2013) estimates should not be trusted. Employing a minimum misfit
suggested a lower limit of detectability of M ∼ 5 at regional criteria into the joint DecisionModule would help throw out
distances, with dynamic motions visible at M 4.6 in the near unreliable source estimates from geodetic algorithms.
field; however, the seismic algorithms will effectively work to Network hardening, in the areas of latency and telemetry
M < 3. Understanding how to weight magnitude estimates robustness, is very important for early warning stability. The
from G-FAST within the joint DecisionModule will be crucial simulations for dropouts showed the greatest variability for the
going forward. The ShakeAlert DecisionModule takes the out- PGD depth and magnitude estimates and indicate that improv-
put of all the seismic algorithms and determines a best magni- ing telemetry links provides a high return on investment. Laten-
tude, location, and OT; modifying this to handle PGD cies, in the near field, are vital for reducing the no-warning-zone,
magnitudes/depths and finite-fault information is complicated both from seismic and geodetic algorithms, and can be reduced
and an area of active research, which we discuss in some of the by improving telemetry, reducing packet sizes, or increasing sta-
possibilities below. Grapenthin et al. (2014b) looked at after- tion density.

940 Seismological Research Letters Volume 87, Number 4 July/August 2016

 5 (i.e., MMI between VI and VIII); however, how this would
be perceived by the general public is a complex issue, and lack
of decisiveness will be considered a system weakness.
4 Current Station Configuration
Δ VR (Main − Auxiliary) (%)

CONCLUSIONS
3 We have shown the utility of geodetic EEW for a moderate-
sized deep earthquake under Puget Sound. G-FAST is capable
of providing timely, stable, and robust estimates of magnitude,
2 depth, and source parameters under real-world conditions. Es-
timates of magnitude and depth from PGD scaling are quick
and accurate; solutions for this simulation of the Nisqually
Old Station Configuration earthquake are available by 22 s after OT, considering latency,
1
and fully stable by 30 s, converging to the correct magnitude
and depth. The CMT results are available by 43 s after OT
when considering latency and stabilize after ∼50 s. The CMT
0
0 250 500 750 1000
nodal plane results are within 10° in both strike and dip from
Iteration the input model and the main fault-plane solution from the
USGS and within several kilometers in depth. The final-mag-
▴ Figure 9. Difference in VR between inversion on main and aux- nitude estimate from the CMT is 0.3 magnitude units larger
iliary fault planes for 1000 noise simulations with 1 m of dip slip on than postprocessed results. The slip inversion was not capable
the steeply dipping fault plane (main). The inversion results using of ascertaining the correct fault plane for this deep event; how-
the current (2015) station configuration are shown in black and ever, with the current network configuration, we expect im-
the old station configuration (2001; strong-motion stations) used in provements in this capability. Previous studies have shown
this study is shown in gray. the ability for geodetic EEW to work for all earthquakes of soci-
etal relevance at regional distances, demonstrating that inte-
grating high-rate geodetic observations into EEW systems
worldwide is crucial for a complete and robust system.
The most pressing issue going forward will be the seamless
combination of magnitude and source estimates from the seis-
mic algorithms (ElarmS, OnSite, VS, and FinDer [Böse et al., DATA AND RESOURCES
2012]) and the geodetic algorithms (G-FAST, G-larmS, and
BEFORES) and how to optimize the DecisionModule within Synthetic seismograms used in this study are available via re-
ShakeAlert to use the most appropriate earthquake model. For quest to the corresponding author B. W. C. Seismic data used
example, the magnitude estimates from the seismic algorithms can be obtained from the Incorporated Research Institutions
should be trusted with high confidence for lower-magnitude for Seismology (IRIS) Data Management Center (www.iris.
events, but the geodetic modules should be trusted more for edu). Plots were created using the Generic Mapping Tools
higher magnitudes. v.4.5.9 (www.soest.hawaii.edu/gmt; Wessel and Smith, 1998).
One must be careful, however, not to inadvertently bias The details of the following resources can be found at modified
the source characterization; if magnitude saturation occurs Mercalli intensity (MMI) VI–VII throughout the Puget Low-
within the seismic algorithms, they will report a smaller mag- lands (http://earthquake.usgs.gov/earthquakes/shakemap/pn/
nitude and hence a higher weight than it should deserve in that shake/01022818540, JavaScript Object Notation JSON proto-
instance. In addition, the use of finite faults versus point col (http://www.json.org ), and Python ObsPy package (http://
sources will drastically change the ground-motion prediction www.obspy.org ). All the websites were last accessed January
within the UserDisplay. The UserDisplay takes the magnitude, 2015.
location, and timing information from the DecisionModule
and computes the predicted ground motion and arrival time
at the user’s location. One example of a potential modification ACKNOWLEDGMENTS
to the DecisionModule is to operate on each methodology
independently and then run a weighted average of the ground- We would like to thank John Langbein and Diego Melgar for
motion prediction from each method. This will still introduce thoughtful comments that improved the article. We also thank
biases to the final solution, so a detailed study of these effects two anonymous reviewers and the Editor Zhigang Peng for
would have to be performed to understand the potential im- helpful reviews. This work has been funded by the Gordon
pacts on EEW. One could also simply report a range of poten- and Betty Moore Foundation Grant Number 663450 to Uni-
tial motions at a given location from all possible algorithms versity of Washington.

Seismological Research Letters Volume 87, Number 4 July/August 2016 941

REFERENCES Genrich, J. F., and Y. Bock (2006). Instantaneous geodetic positioning
with 10–50 Hz GPS measurements: Noise characteristics and
Allen, R. M., and A. Ziv (2011). Application of real-time GPS to earth- implications for monitoring networks, J. Geophys. Res. 111,
quake early warning, Geophys. Res. Lett. 38, L16310, doi: 10.1029/ no. B03403, doi: 10.1029/2005JB003617.
2011GL047947. Given, D. D., E. S. Cochran, T. Heaton, E. Hauksson, R. Allen, P.
Altamimi, Z., L. Métivier, and X. Collilieux (2012). ITRF2008 plate mo- Hellweg, J. Vidale, and P. Bodin (2014). Technical implementation
tion model, J. Geophys. Res. 117, no. B07402, doi: 10.1029/ plan for the ShakeAlert production system: An earthquake early
2011JB008930. warning system for the west coast of the United States, U.S. Geol.
Blewitt, G., C. Kreemer, W. C. Hammond, H.-P. Plag, S. Stein, and E. Surv. Open-File Rept. 2014-1097, 25 pp., doi: 10.3133/ofr20141097.
Okal (2006). Rapid determination of earthquake magnitude using Grapenthin, R., I. A. Johanson, and R. M. Allen (2014a). Operational
GPS for tsunami warning systems, Geophys. Res. Lett. 33, L11309, real-time GPS-enhanced earthquake early warning, J. Geophys. Res.
doi: 10.1029/2006GL026145. 119, no. 10, 7944–7965, doi: 10.1002/2014JB011400.
Bock, Y., D. Melgar, and B. W. Crowell (2011). Real-time strong-motion Grapenthin, R., I. Johanson, and R. M. Allen (2014b). The 2014 M w 6.0
broadband displacements from collocated GPS and accelerometers, Napa earthquake, California: Observations from real-time GPS-en-
Bull. Seismol. Soc. Am. 101, no. 6, 2904–2925, doi: 10.1785/ hanced earthquake early warning, Geophys. Res. Lett. 41, no. 23,
0120110007. 8269–8276, doi: 10.1002/2014GL061923.
Böse, M., E. Hauksson, K. Solanki, H. Kanamori, Y.-M. Wu, and T. H. Hartog, J. R., V. C. Kress, S. D. Malone, P. Bodin, J. E. Vidale, and B. W.
Heaton (2009). A new trigger criterion for improved real-time Crowell (2016). Earthquake early warning: ShakeAlert in the
performance of onsite earthquake early warning in southern Cal- Pacific Northwest, Bull. Seismol. Soc. Am. (accepted).
ifornia, Bull. Seismol. Soc. Am. 99, no. 2A, 897–905, doi: 10.1785/ Hashima, A., Y. Takada, Y. Fukahata, and M. Matsu’ura (2008). General
0120080034. expressions for internal deformation due to a moment tensor in an
Böse, M., T. H. Heaton, and E. Hauksson (2012). Real-time finite elastic/viscoelastic multilayered half-space, Geophys. J. Int. 175,
fault rupture detector (FinDer) for large earthquakes, Geophys. J. no. 3, 992–1012, doi: 10.1111/j.1365-246X.2008.03837.x.
Int. 191, no. 2, 803–812, doi: 10.1111/j.1365-246X.2012.05657.x. Ichinose, G. A., H. K. Thio, and P. G. Somerville (2004). Rupture process
Böse, M., T. Heaton, and K. Hudnut (2013). Combining real-time and near-source shaking of the 1965 Seattle–Tacoma and 2001
seismic and GPS data for earthquake early warning (invited), Nisqually, intraslab earthquakes, Geophys. Res. Lett. 31, L10604,
American Geophysical Union, Fall Meeting 2013, abstract number doi: 10.1029/2004GL019668.
G51B-05. Kao, H., K. Wang, R.-Y. Chen, I. Wada, J. He, and S. D. Malone (2008).
Brown, H. M., R. M. Allen, M. Hellweg, O. Khainovski, D. Neuhauser, Identifying the rupture plane of the 2001 Nisqually, Washington,
and A. Souf (2011). Development of the ElarmS methodology for earthquake, Bull. Seismol. Soc. Am. 98, no. 3, 1546–1558, doi:
earthquake early warning: Realtime application in California and 10.1785/0120070160.
offline testing in Japan, Soil Dyn. Earthq. Eng. 31, no. 2, 188– Kuyuk, H. S., R. M. Allen, H. Brown, M. Hellweg, I. Henson, and D.
200, doi: 10.1016/j.soildyn.2010.03.008. Neuhauser (2014). Designing a network-based earthquake early
Bustin, A., R. D. Hyndman, A. Lambert, J. Ristau, J. He, H. Dragert, and warning algorithm for California: ElarmS-2, Bull. Seismol. Soc.
M. Van der Kooij (2004). Fault parameters of the Nisqually earth- Am. 104, no. 1, 162–173, doi: 10.1785/0120130146.
quake determined from moment tensor solutions and the surface Melgar, D., Y. Bock, and B. W. Crowell (2012). Real-time centroid
deformation from GPS and InSAR, Bull. Seismol. Soc. Am. 94, moment tensor determination for large earthquakes from local and
no. 2, 363–376, doi: 10.1785/0120030073. regional displacement records, Geophys. J. Int. 188, no. 2, 703–718,
Colombelli, S., R. M. Allen, and A. Zollo (2013). Application of real-time doi: 10.1111/j.1365-246X.2011.05297.x.
GPS to earthquake early warning in subduction and strike-slip envi- Melgar, D., Y. Bock, D. Sanchez, and B. W. Crowell (2013). On robust
ronments, J. Geophys. Res. 118, no. 7, 3448–3461, doi: 10.1002/ and reliable automated baseline corrections for strong motion seismol-
jgrb.50242. ogy, J. Geophys. Res. 118, no. 3, 1177–1187, doi: 10.1002/jgrb.50135.
Crowell, B. W., Y. Bock, and D. Melgar (2012). Real-time inversion of Melgar, D., B. W. Crowell, Y. Bock, and J. S. Haase (2013). Rapid mod-
GPS data for finite fault modeling and rapid hazard assessment, eling of the 2011 M w 9.0 Tohoku-Oki earthquake with seismogeod-
Geophys. Res. Lett. 39, L09305, doi: 10.1029/2012GL051318. esy, Geophys. Res. Lett. 40, no. 12, 2963–2968, doi: 10.1002/
Crowell, B. W., Y. Bock, and M. Squibb (2009). Demonstration of earth- grl.50590.
quake early warning using total displacement waveforms from real Melgar, D., B. W. Crowell, J. Geng, R. M. Allen, Y. Bock, S. Riquelme, E.
time GPS networks, Seismol. Res. Lett. 80, no. 5, 772–782, doi: Hill, M. Protti, and A. Ganas (2015). Earthquake magnitude calcu-
10.1785/gssrl.80.5.772. lations without saturation from the scaling of peak ground displace-
Crowell, B. W., D. Melgar, Y. Bock, J. S. Haase, and J. Geng (2013). ment, Geophys. Res. Lett. 42, no. 13, 5197–5205, doi: 10.1002/
Earthquake magnitude scaling using seismogeodetic data, 2015GL064278.
Geophys. Res. Lett. 40, no. 23, 6089–6094, doi: 10.1002/ Minson, S. E., J. R. Murray, J. O. Langbein, and J. S. Gomberg (2014).
2013GL058391. Real-time inversions for finite fault slip models and rupture geom-
Cua, G., and T. Heaton (2007). The virtual seismologist (VS) method: etry based on high-rate GPS data, J. Geophys. Res. 119, no. 4, 3201–
A Bayesian approach to earthquake early warning, in Earthquake 3231, doi: 10.1002/2013JB010622.
Early Warning Systems, P. Gasparini, G. Manfredi, and J. Zschau Ohta, Y., T. Kobayashi, H. Tsushima, S. Miura, R. Hino, T. Takasu, H.
(Editors), Springer, Berlin, Germany, 97–130, doi: 10.1007/978- Fujimoto, T. Iinuma, K. Tachibana, T. Demachi, et al. (2012).
3-540-72241-0_7. Quasi real-time fault model estimation for near-field tsunami fore-
Dreger, D., and A. Kaverina (2000). Seismic remote sensing for the earth- casting based on RTK-GPS analysis: Application to the 2011 To-
quake source process and near-source strong shaking: A case study hoku-Oki earthquake (M w 9.0), J. Geophys. Res. 117, no. B02311,
of the October 16, 1999 Hector Mine earthquake, Geophys. Res. doi: 10.1029/2011JB008750.
Lett. 27, no. 13, 1941–1944. Okada, Y. (1985). Surface deformation to shear and tensile faults in a
Geng, J., Y. Bock, D. Melgar, B. W. Crowell, and J. S. Haase (2013). A half-space, Bull. Seismol. Soc. Am. 75, 1135–1154.
new seismogeodetic approach applied to GPS and accelerometer O'Toole, T. B., A. P. Valentine, and J. H. Woodhouse (2013). Earthquake
observations of the 2012 Brawley seismic swarm: Implications for source parameters from GPS-measured static displacements with
earthquake early warning, Geochem. Geophys. Geosyst. 14, no. 7, potential for real-time application, Geophys. Res. Lett. 40, no. 1,
2124–2142, doi: 10.1002/ggge.20144. 60–65, doi: 10.1029/2012GL054209.

942 Seismological Research Letters Volume 87, Number 4 July/August 2016

 Wessel, P., and W. H. F. Smith (1998). New, improved version of generic Earthquake Science Center
mapping tools released, Eos Trans. AGU 79, no. 47, 579, doi: University of Washington
10.1029/98EO00426.
Wright, T. J., N. Houlié, M. Hildyard, and T. Iwabuchi (2012). Real- Johnson Hall Room-070, Box 351310
time, reliable magnitudes for large earthquakes from 1 Hz GPS pre- 4000 15th Avenue NE
cise point positioning: The 2011 Tohoku-Oki (Japan) earthquake, Seattle, Washington 98195-1310 U.S.A.
Geophys. Res. Lett. 39, L12302, doi: 10.1029/2012GL051894.
Zhu, L., and L. A. Rivera (2002). A note on the dynamic and static displace- Timothy I. Melbourne
ments from a point source in multi-layered media, Geophys. J. Int. 148,
no. 3, 619–627, doi: 10.1046/j.1365-246X.2002.01610.x. Marcelo Santillan
Zumberge, J. F., M. B. Heflin, D. C. Jefferson, M. M. Watkins, and F. H. Department of Geological Sciences
Webb (1997). Precise point positioning for the efficient and robust Central Washington University
analysis of GPS data from large networks, J. Geophys. Res. 102, 400 E. University Way, MS 7418
no. B3, 5005–5017, doi: 10.1029/96JB03860. Ellensburg, Washington 98926 U.S.A.
Brendan W. Crowell Sarah E. Minson
David A. Schmidt U.S. Geological Survey
Paul Bodin Earthquake Science Center
John E. Vidale 345 Middlefield Road, MS 977
J. Renate Hartog Menlo Park, California 94025-3591 U.S.A.
Victor C. Kress
Department of Earth and Space Sciences Dylan G. Jamison
University of Washington Department of Earth and Environmental Sciences
Johnson Hall Room-070, Box 351310 University of Waterloo
4000 15th Avenue NE 200 University Avenue West
Seattle, Washington 98195-1310 U.S.A. Waterloo, Ontario
crowellb@uw.edu Canada N2L 3G1
Joan Gomberg Published Online 8 June 2016
U.S. Geological Survey

Seismological Research Letters Volume 87, Number 4 July/August 2016 943

View publication stats