# Earthquake nucleation and fault slip complexity in the lower crust of central Alaska

Carl Tape, Stephen Holtkamp, Vipul Silwal, Jessica Hawthorne, Yoshihiro

Kaneko, Jean Paul Ampuero, Chen Ji, Natalia A Ruppert, Kyle Smith,

Michael West

► To cite this version:

Carl Tape, Stephen Holtkamp, Vipul Silwal, Jessica Hawthorne, Yoshihiro Kaneko, et al.. Earthquake nucleation and fault slip complexity in the lower crust of central Alaska. Nature Geoscience, 2018, 11 (7), pp.536-541. 10.1038/s41561-018-0144-2. hal-02019296

HAL Id: hal-02019296

https://hal.science/hal-02019296v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Earthquake nucleation and fault slip complexity in the
lower crust of central Alaska

Carl Tape<fnref n="1" />, Stephen Holtkamp<fnref n="1" />, Vipul Silwal<fnref n="1" />, Jessica Hawthorne<fnref n="2" />, Yoshihiro Kaneko<fnref n="3" />, Jean Paul Ampuero<fnref n="4" />, Chen Ji<fnref n="5" />, Natalia Ruppert<fnref n="1" />, Kyle Smith<fnref n="1" />, & Michael E. West<fnref n="1" />

Earthquakes start under conditions that are largely unknown. In laboratory analog experi-
ments and continuum models, earthquakes transition from slow-slipping, growing nucleation
to fast-slipping rupture. In nature, earthquakes generally start abruptly, with no evidence
for a nucleation process. Here we report evidence from a strike-slip fault zone in central
Alaska of extended earthquake nucleation and of very-low-frequency earthquakes (VLFE),
a phenomenon previously reported only in subduction zone environments. In 2016 a VLFE
transitioned into an earthquake of magnitude 3.7 and was preceded by a 12-hour-long ac-
celerating foreshock sequence. Benefiting from 12 seismic stations deployed within 30 km of
the epicenter, we identify coincident radiation of distinct high-frequency and low-frequency
waves during 22 s of nucleation. The power-law temporal growth of the nucleation signal
is quantitatively predicted by a model in which high-frequency waves are radiated from the

<footnote n="1">¹Geophysical Institute, University of Alaska Fairbanks, Fairbanks, AK, 99775, USA</footnote>

<footnote n="2">²Department of Earth Sciences, University of Oxford, Oxford, England OX1 3AN, UK</footnote>

<footnote n="3">³GNS Science, Lower Hutt, 5040, New Zealand</footnote>

<footnote n="4">⁴Seismological Laboratory, California Institute of Technology, Pasadena, CA, 91125, USA</footnote>

<footnote n="5">⁵Department of Earth Science, University of California, Santa Barbara, CA, 93106, USA</footnote>
vicinity of an expanding slow slip front. The observations reveal the continuity and complexity of slip processes near the bottom of the seismogenic zone of a strike-slip fault system in central Alaska.

Large faults accommodate relative plate motion through a wide variety of processes ranging from earthquakes and aseismic transient slip events to stable sliding¹. Earthquakes occur within subducting slabs or are confined to the shallowest, brittle regimes of Earth's crust, usually no deeper than 15–20 km for crustal or strike-slip fault systems or 30–50 km for subduction zones. Aseismic transient slip events tend to cluster near the frictionally unstable-to-stable transition¹, while stable sliding occurs deeper, where the fault zone is frictionally stable or ductile at plate convergence strain rates.

Earthquake occurrence has traditionally been viewed as a stochastic process, but observations of foreshock sequences² have suggested that large earthquakes might be preceded by an extended nucleation process³⁻⁶ which, if recognized, could be used as a predictive measure³,⁷,⁸. Unfortunately, foreshock sequences have only been recognized as such after the fact and are not observed in many earthquakes.

The discovery of aseismic transient events at plate margins worldwide¹ has led to many theories and observations attempting to relate slow earthquake phenomena to large earthquake rupture⁴,⁹⁻¹¹. Slow earthquake observations, which include geodetically observable slow-slip events¹ and seismically observable slow earthquakes such as low-frequency earthquakes (LFEs) and very-low-frequency earthquakes (VLFEs)¹², have provided a physical basis for some foreshock sequences⁴.
However, geodetic and seismic observations generally provide a limited view of precursory slow slip. For example, though it is recognized that aseismic slip preceded the 2011 Tohoku-oki earthquake<sup>13</sup>, the inferred slip is limited to the weeks to days prior to the mainshock rupture. This leads to an observational gap in the nucleation process: What, if anything, is occurring in the minutes to seconds before mainshock rupture?

Significant efforts to identify a nucleation process have been attempted in theoretical studies<sup>14</sup>, in numerical modeling<sup>15-18</sup>, and in laboratory experiments<sup>19-23</sup>, but they have proven challenging in observational seismology<sup>4,24,25</sup>. Laboratory experiments<sup>23</sup> show that earthquake nucleation exhibits an early stage, characterized by low (although increasing) rupture velocities until rupture approaches a critical size, and a later stage, where rupture velocity rapidly increases. Precursory aseismic transients, observed geodetically<sup>26,27</sup> and interpreted in many cases as drivers of earthquake foreshock swarms<sup>13</sup>, may represent an early stage of earthquake nucleation. Attempts have been made to seismically observe the later stage of nucleation<sup>4,24,28</sup>. Near-source borehole strain-meter data from California and Japan have detected slow earthquake transients<sup>27,29</sup> but have not revealed precursory signals for earthquakes, suggesting that a nucleation process is not common (or observable) for these fault zones. The continuum of nucleation processes has not been observed in one setting, for example, by a slow earthquake—or a cascade of smaller events—directly transitioning into a normal earthquake.

In 2012, Tape et al. [28] identified a signal of earthquake nucleation that lasted 24 seconds before transitioning into a magnitude 3.8 earthquake. The earthquake source mechanism was con-
sistent with left-lateral strike-slip faulting within the Minto Flats fault zone, a 180-km-long fault zone with two main faults, one of which produced a magnitude 6.0 earthquake in 1995 (Fig. 1)³⁰,³¹. Interpretation of the 2012 event as nucleation+earthquake was complicated by the simultaneous arrival of surface waves from a $M_w$ 8.6 earthquake in the Indian Ocean³². It seemed clear that the $M_w$ 8.6 surface waves triggered the nucleation, but it was not clear that the nucleation and subsequent earthquake would have occurred in the absence of the surface waves.

An installation of 13 seismometers in the Minto Flats fault zone in 2015 led to the discovery of VLFEs, reported next, one of which transitioned abruptly into a typical earthquake. We describe two unusual events that occurred following the installation.

## Very low frequency earthquakes in a strike-slip fault

The first event occurred on 2015-09-12 and was identified on dozens of stations in central Alaska. The waveforms are characterized by an emergent signal and do not exhibit the distinctive P and S waves of a typical earthquake (Fig. 2a). However, a weak but identifiable impulse, likely a P-phase of an initiating earthquake, occurred at the beginning of the event. Analysts estimated a magnitude of 2.6 based on the P amplitudes. By filtering the seismogram at high frequencies and then plotting the envelope of the log-scaled amplitudes, we found additional complexity: the event was preceded by a smaller, similar event 70 seconds earlier.

The long-period (>10 s) filtered seismograms of the 2015 event reveal a VLFE (Fig. 2a), the first observed in an intraplate setting and the first associated with a strike-slip fault system. A
remarkably simple set of waveforms is visible on stations out to distances of 1000 km (Fig. S1).

The waveforms provide high-quality estimates for the magnitude ($M_w$ 3.8), centroid epicentral location, depth (21 km; Fig. S2), and mechanism of the event (Fig. S3). The source duration for a typical earthquake of this magnitude is about 0.5 s, but the 2015 event lasted 10 s (Figs. S4-S5).

Compared to a normal earthquake with similar amplitude at high frequencies and similar depth, the 2015 event is enhanced at low frequencies (Fig. S6), in a similar fashion to VLFEs in subduction settings$^{33,34}$.

High-frequency and low-frequency waves from VLFEs have been interpreted to represent two facets of the same process, with two possible interpretations, both involving numerous small events radiating high-frequency waves. In one interpretation, there is a slow-slipping process emanating low-frequency waves, with numerous small events at the margins or inside of the slow-slipping pulse$^{35,36}$. An alternative interpretation is that the waves from the small events interfere in a manner that produces the appearance of a simple low-frequency signal$^{37,38}$. Within subduction settings, VLFEs are associated with LFEs and tectonic tremor, which are best observed in the ~2–8 Hz range$^{1,12}$. The high-frequency components associated with the 2015 VLFE are not band-limited, and the spectra look similar to those of normal earthquakes, except for the characteristic enhancement at long periods (Fig. S6).

It is hard to say when the 2015 event started, since the P onset appears within the coda of an earlier, much smaller event (Fig. 2a). Differential travel times and high inter-station phase coherence$^{39}$ between both signals imply that the earlier event originated from the same hypocenter
and with the same mechanism as the VLFE. The similarity in high-frequency envelopes and in
spectra suggests that the earlier event was also a VLFE but that its long-period signal was below
noise levels. It is plausible that the 2015 event initiated as a small VLFE lasting ~10 s, and then,
after about 60 s, transitioned into a larger VLFE at 21 km depth that slipped for about 10 seconds
in a left-lateral sense on the Minto Flats fault zone. The P phase at the onset of the larger VLFE
suggests that fast-slipping sources are intertwined with the slow-slip progression.

Finding the 2015 VLFE prompted us to look, though not exhaustively, for other VLFEs in
previous years using regional long-period data. We identified a second, smaller ($M_w$ 3.5) VLFE in
2013 (Fig. 1b) with a source duration of 12 s. As with the 2015 event, the source mechanism is
consistent with left-lateral faulting at 23 km depth (Fig. 1a). The occurrence of the 2013 and 2015
VLFEs, separated by 50 km, suggests that a similar deep slow-slipping process may be possible
throughout the Minto Flats fault zone.

VLFE precursor to an earthquake

We now turn to the second event, a magnitude 3.7 earthquake that occurred on 2016-01-14 at a
depth of 17 km with an epicenter near the center of our array of stations (Fig. 1a). The earthquake
had clear P waves, S waves, and sufficient low-frequency waves to allow for a source mechanism
to be determined (Fig. S3b). The mechanism is consistent with left-lateral faulting. Given its
magnitude, the source dimension is between 0.5 km and 2.5 km, based on an assumed stress drop
between 0.2 and 20 MPa [⁴⁰].
The envelopes of the high-frequency seismograms of the 2016 event reveal a growing nucle-
ation signal lasting about 22 s (Fig. 2b), reminiscent of the 2012 triggered earthquake (Fig. 1b). In
this case, however, there was no remote earthquake, and we have many stations much closer to the
epicenter: 12 within 30 km. We can therefore examine subtle signals associated with nucleation.

A close examination of the low-frequency (<0.1 Hz) seismograms reveals two parts: (1) a low-amplitude signal occurring during the high-frequency foreshock signal and (2) a high-amplitude signal occurring during the main earthquake. The novel signal is the first one, plotted at the top of Figure 2b, and identified on 11 of the 12 closest stations. The low-frequency foreshock signal clearly arrives at the closest stations before the origin time of the mainshock. Based on the similarity to the 2015 VLFE (emergent high frequency waveforms and coincident low frequency signal), we classify this low-frequency foreshock as a VLFE.

Although we do not capture enough of the low-frequency foreshock waveforms to directly
invert for a source mechanism, we are able to compare three sets of polarity measurements for
the 2016 event (Tab. S3): the polarity of the low-frequency foreshock signal, the polarity of the
P wave for the Mw 3.7 earthquake, and the predicted polarity from the mechanism inferred from
the long-period waves. The results provide strong evidence for the same left-lateral slip process
occurring during the nucleation stage and rupture of the earthquake.

The exceptional station coverage for the 2016 event allows us to demonstrate that the two
distinct signals of the VLFE and the earthquake (Fig. 2b) originate from the same location and
source mechanism, and are therefore likely to be part of the same continuous process. We apply
a new technique³⁹ to show that the foreshock signal is no more than 1 km from the mainshock signal. The technique compares 4-second time windows of high-frequency waveforms from the mainshock with waveforms from the foreshock signal. We find that the inter-station coherence is high between 18 and 1.5 s before the mainshock (Fig. 3). A grid search of possible foreshock locations gives the highest coherence when the foreshocks are assumed to be 0–1 km deeper than the mainshock centroid, suggesting that the VLFE occurred in the transition zone to more frictionally stable rocks.

While the high-frequency foreshock signal and precursory VLFE occurred tens of seconds prior to the 2016 mainshock, the nucleation process may have started hours beforehand. A catalog of earthquakes within 5 days and ~20 km of the mainshock (Fig. S7) shows that an accelerating occurrence of foreshocks began about 12 hours prior to the mainshock. Accelerating foreshocks have been observed prior to large earthquake and interpreted as aseismic fault slip⁴, ⁶, ¹³.

## Interpretations of fault slip processes

Our observations from a 60-km section of the Minto Flats fault zone (Fig. 1) provide evidence for long-duration earthquakes (VLFEs), ‘normal’ fast-slip earthquakes, and events that display the transition between the two. The events in the fault zone, at depths between 16 km and 24 km, provide a window into slip processes near the base of the seismogenic zone.

Our observations allow us to identify and discuss two stages for long-duration events in the fault zone: (1) nucleation, manifested by ~20 seconds of simultaneous high-frequency foreshocks
and a low-frequency signal of a VLFE; and (2) an earthquake initiated by that VLFE. Five events in this study are summarized in Table S2. Based on observations of the 2016 event, we now favor an interpretation for the 2012 event²⁸ as a VLFE—evidenced by the 24 s foreshock signal in Figure 1b—that transitioned into a M<sub>w</sub> 3.8 earthquake. The VLFEs of 2013 and 2015 exhibited the nucleation stage but did not lead to earthquakes. The events of 2012 and 2016 spanned both stages. Our observations suggest that the high-frequency nucleation signal will coincide with VLFEs, which in turn may or may not transition into earthquakes.

We consider different mechanisms for explaining the two stages of nucleation and earthquake, as shown in the cartoon in Figure 4, which is based on high-quality seismic observations as well as on modeling efforts (see Methods). For each of the two stages, we consider two possibilities, leading to four scenarios (*a* → *b*, *a* → *d*, *c* → *b*, *c* → *d*) for the 2016 event. For the M<sub>w</sub> 3.7 earthquake of the 2016 event, its dynamic rupture could transition directly from a slow slip event (Fig. 4b) or it could be triggered from a local event, without slow slip (Fig. 4d). For the nucleation stage, the combined high-frequency and low-frequency signals could either be a slow-slipping event with associated high-frequency bursts (Fig. 4a), analogous to aseismic slow slip associated with tremor<sup>41</sup>, or they could be a cascading sequence of small events<sup>24,37</sup>, without underlying slow slip, that collectively generate the appearance of a long-period signal<sup>38</sup> (Fig. 4c).

Smaller events are present within the longer-duration events in our study, as evidenced from phase coherence analysis, source time functions, and analysis of spectra. Phase coherence modeling of observations (Fig. 3) for the 2016 event showed that the high-frequency foreshocks orig-
inated within 1 km of the $M_w$ 3.7 earthquake, similar to the case of the 2012 event³⁹. High-frequency ($f \le 4.5$ Hz) modeling of waveforms produced a source time function with multiple peaks, representing distinct subevents. Phase coherence within the 2015 VLFE (and also to the preceding, smaller VLFE) revealed the presence of smaller events. The source time function estimated for the 2015 VLFE (Fig. S4a) revealed 8–10 subevents within the 10-s VLFE duration. This complex source time function, along with an analysis of the displacement spectra (Section S1), imply that it is possible that most of the VLFE slip accumulates in ~10 closely spaced $M_w$ 3.2 earthquakes. However, it is exceedingly uncommon to observe 10 $M_w$ 3.2 earthquakes in 10 seconds without some underlying mechanism, such as pore fluid diffusion⁴², rock damage⁴³, or aseismic slip¹⁴⁻²³.

Aseismic slip is an appealing mechanism—which could accommodate much of the moment in the 2015 and 2016 events—because the evolution of seismic amplitude prior to the 2016 earthquake matches the growth predicted by frictional models of earthquake nucleation. For the 2016 event the high frequency foreshock signals lasted 22 s and grew according to $t^k$, where $t$ is time and where the exponent $k$ is estimated from different stations to be between 2 and 4 (Fig. S8). Similar exponents are found for the spatial growth of expanding earthquake nucleation fronts in recent numerical modeling⁴⁴ and laboratory observations²³ (Fig. 5). A plausible interpretation is that the growing seismic amplitudes reflect a number of small foreshocks being triggered as the slow rupture propagates outward (see Methods).

Aseismic slip preceding the 2016 VLFE and earthquake is supported by observations of the
~12-hour accelerating foreshock sequence (Fig. S7). Globally, large earthquakes have been preceded by accelerating foreshock sequences and geodetically observed slow slip<sup>5,9,26</sup>. For the 2016 event it seems plausible that aseismic slow slip transitioned into a seismically observable VLFE with associated high-frequency foreshocks (Fig. 4a). Unlike the 2016 event, the 2015 VLFE did not exhibit an accelerating foreshock sequence. Thus other VLFEs may start like most earthquakes: without any warning. The 2012 triggered VLFE was preceded by a small swarm of ~10 earthquakes, indicating that aseismic processes may have already been underway when the large surface waves triggered the VLFE, and the 2013 VLFE was preceded by two possible foreshocks.

Recent laboratory experiments and numerical models<sup>45-47</sup> imply that the slow-slipping (Fig. 4a) and cascading (Fig. 4c) nucleation processes may be concurrent, with both playing a role. In this regard, the two scenarios could be thought of as end members. For two events presented here, the nucleation process was manifest in novel observations lasting ~20 s (24 s in 2012; 22 s in 2016) before transitioning into an earthquake.

The Minto Flats fault zone in central Alaska, capable of $M > 6$ earthquakes, is a natural laboratory for studying earthquake nucleation and the continuum of fault rupture processes that can occur on large-scale faults. With multiple observations distributed along the fault zone, the MFFZ appears to have the right specialized conditions to generate extended nucleation sequences, unlike other well instrumented faults such as the San Andreas Fault where earthquake precursors have been searched for unsuccessfully. From our modeling, the size of the nucleation zone is proportional to the critical nucleation length, which depends on the effective normal stress
and temperature-controlled frictional properties of fault-zone rocks<sup>44,48</sup>. Critical nucleation lengths would be larger near the base of the seismogenic zone due to high fluid pressure or velocity-neutral frictional properties of the fault-zone rocks. Therefore the transitional region is prone to generating nucleation sequences with large aseismic moment and numerous accompanying foreshocks. Foreshocks (localized rapid slip) during slow nucleation could be promoted by heterogeneity deep within the fault zone. In a view that has emerged from a decade of observations of tectonic tremor during slow slip¹ and deep high-frequency radiation during large earthquakes<sup>49</sup>, fault zone materials in the deep transition zone are thought to be highly heterogeneous, with intermingled seismic and aseismic behavior.

We have observed a nucleation process that has previously been limited to laboratory experiments<sup>20,23</sup>, where slip processes are examined within cm-scale faults creating earthquakes with magnitudes less than −4 [<sup>22,23</sup>]. The events identified in our study, with magnitudes 3.4–3.9, provide an observational bridge from laboratory to real world environments.

## Methods

Methods, including statements of data availability and any associated accession codes and references, are available in the online version of this paper.

**Seismogram processing.** Instrument responses were deconvolved from all seismograms using Seismic Analysis Code. The details of deconvolution and subsequent filtering and processing are important for our analysis, because we are discussing precursory low-frequency signals. To
avoid any potential for precursory ringing due to filtering, we employ only one-pass causal filters
throughout instrument deconvolution and subsequent processing. As a result, the sign of the phase
and onset timing of arriving low-frequency signals are preserved (e.g., arriving low-frequency
signals will not be arriving early due to filtering).

We examine the onset of the precursory low-frequency signal, which we identify as a VLFE,
to the 2016 earthquake (Fig. 2, Figs. S10–S12). Using causal filters and in the presence of noise,
we expect the picked arrivals to be later than the actual arrivals of the low-frequency waves. With
Figure S9, we attempt to determine the amount of time lag expected in the presence of noise. We
create an artificial signal (an upward pointing triangle function beginning at t = 0 and with a
width of 20 seconds), add Gaussian noise, apply a causal bandpass filter, and pick the onset of the
resulting low-frequency signal above the noise level. From 10% to 40% SNR noise levels, the pick
times range from 8.3 to 11.2 seconds, suggesting that the low-frequency onset time we picked is
likely delayed by ~8–12 seconds from the actual onset time of the signal.

**Source mechanism inversions and estimation of source duration.** We estimated source mechanisms using the method of ref ⁵⁰, which uses five time windows of waveforms for each station: P on the vertical component, P on the radial component, Rayleigh on the vertical component, Rayleigh on the radial component, and Love on the transverse component. Different bandpass filters are applied for body waves and surface waves in order to maximize the signal-to-noise ratios.

Figure S3 shows fits between observed seismograms and synthetic seismograms computed
using the source mechanisms indicated by the beachball focal mechanisms. The estimated depth
for each event is shown in Figure S2.

Careful examination of signal-to-noise ratios allowed us to isolate some of the unusual events and processes in this study. Figure S1 shows a record section for the 2015 VLFE, filtered showing clear long-period signals (20–50 s) for stations out to 1000 km from the epicenter. Figure S6 shows signal and noise spectra for the 2015 VLFE and a comparison earthquake. The spectra show how the VLFE is characterized by amplified seismic energy at periods above 10 s.

We estimated the source duration of the events using two approaches. In the first approach, we allowed the source duration to be varied within the source mechanism inversion. For the 2015 VLFE, the best-fitting synthetics were generated using a source duration that was much longer than expected for an earthquake with the same magnitude ($M_w$ 3.8). In the second approach, we use each P and SH wave for each station to estimate a source-time function (Fig. S4) using Landwebber deconvolution⁵¹. Stations were selected based on waveform fits, using a criteria of 50% variance reduction when using the estimated source time function versus a delta function. Both approaches yielded a duration of about 10 s for the 2015 event.

**Envelopes of high-frequency seismograms.** High-frequency seismograms are conveniently analyzed using log-scaled envelopes of the seismograms ²⁸,⁵². The processing steps to calculate the envelopes $E(t)$ can be summarized as

$$E(t) = S_N [\log_{10} |\mathcal{H}(F[v(t)])|], \quad (1)$$

where $v(t)$ is the velocity seismograms with instrument response removed, $F$ represents a bandpass filter over the range 2–8 Hz, $\mathcal{H}$ is a Hilbert transform, and $S$ is a N-point median smoothing filter
(we use $N = 10$).

Examples of envelopes are shown for station MDM for several different events in Figure 1b.

Laboratory experiments and numerical experiments offer insights as to the possible functional relationships that may occur during nucleation. We consider a variety of functions in order to fit the observed log-scaled high-frequency envelopes for different events and different stations.

Our preferred function for the HFF is a log function,

$$ \log_{10} f(t) = k \log_{10} (t - t_0) + b \quad (2) $$

which is equivalent to

$$ f(t) = 10^b (t - t_0)^k. \quad (3) $$

Examples of curves of Eq. 2 are shown for all 34 stations in Figure S16. The values of the exponent $k$ cluster between 2 and 4 for stations within 80 km of the epicenter (Fig. S8).

**Estimating the hypocenter for the high-frequency foreshock signal.** We use a matched field technique based on inter-station phase coherence³⁹. The method is designed to identify co-located signals, even if their sources have complex, long-duration source time functions. We can therefore consider a sequence of closely spaced foreshocks as a single complex source, and search for its location relative to the mainshock.

To prepare the data, we bandpass filter between 0.75 and 20 Hz. We extract 5 seconds of the mainshock signal to use as the template, tapering in the first and last 0.5 seconds. For stations farther than 40 km from the source, these 5 seconds start 0.5 seconds before the P arrival. For
more nearby stations, the template starts 0.5 seconds before the expected S arrival, taken to be
1.3 times the P-wave picks, so that each template includes only one of the two arrivals. We then
take the foreshock data—before the mainshock P arrival—and taper the seismograms to zero after
mainshock P arrival, in order to avoid large apparent noise associated with the mainshock.

We cross-correlate the mainshock template with the entire foreshock signal at each station. If the mainshock and foreshocks are co-located, these cross-correlations should eliminate the Green's functions phases, leaving the relative phases of the source time functions, which should be the same across all stations. So we compute the inter-station phase coherence between the cross-correlations in 4-second intervals, averaging over frequencies in the 1–10 Hz band. The calculated phase coherence is larger than the 95% range expected for noise during much of the precursor (Fig. S19a), suggesting that the mainshock and foreshocks are nearly co-located.

To better estimate the allowable foreshock locations, we consider the phase coherence for a
grid of possible foreshock locations within 3 km of the mainshock. For each proposed foreshock
location, we shift the foreshock seismograms by the difference in travel time between that location
and the mainshock, and then compute the phase coherence as described above. The expected
travel times are calculated using the same 1D model that is used by the Alaska Earthquake Center
to locate earthquakes in this region. Figures 3, S20, and S21 illustrate the phase coherence as a
function of the proposed foreshock locations. The highest phase coherence places the foreshocks
< 1 km to the north and 0–1 km below the average mainshock location.
**Accelerating foreshock swarm prior to the 2016 event.** The 2016 event was preceded by about 12 hours of accelerating foreshock occurrence (Fig. S7). We utilize a network matched filter approach to obtain a catalog of seismicity ±5 days and ±20 km of the mainshock. We set a relatively high threshold of 10 times the median absolute deviation, across 39 channels in the 13-station FLATS network (XV), to ensure false positives are not impacting our conclusions. We find that foreshock seismicity started ~12 hours before the mainshock, with a total of ~20 events before the mainshock. For comparison, we repeat this process for a 2014 M<sub>w</sub> 5.0 earthquake that occurred at the northern end of the Minto Flats fault zone, to show a regular mainshock–aftershock sequence, without detectable foreshocks.

**Nucleation phase in a continuum model of a rate-and-state fault.** We analyze rupture growth in a physics-based models of shear-slip nucleation⁴⁴ in order to relate the observed foreshock signals to the potential nucleation process of the 2016 Mₗ 3.7 earthquake. Kaneko et al.⁴⁴ used a two-dimensional (2D) in-plane elastodynamic model⁵³ in which shear-slip events (i.e., earthquakes) are simulated as a part of spontaneously occurring earthquake sequences on a fault embedded into an infinite elastic medium. The model setup was motivated by the laboratory experiments of ref ²³ (Fig. S22a,b). The fault constitutive response is represented by rate-and-state friction laws with the so-called slip law for the state-variable evolution⁵⁴,⁵⁵. The fault is divided up into three segments consisting of a central rate-weakening patch and the surrounding rate-strengthening segments, and is loaded by a background time-independent stressing rate uniformly applied along the fault (Fig. S22b). The detail of the model set-up is described in ref ⁴⁴.

Using this relatively simple model, ref ⁴⁴ reproduced the characteristics of shear-rupture
nucleation as observed in these laboratory experiments (Fig. S22c). In particular, the evolution
of the rupture front under a range of normal stresses σ closely matches the laboratory results
(Fig. S22d,e). In laboratory experiments and in simulations, three distinct phases of rupture evo-
lution were observed regardless of σ: quasi-static propagation characterized by the first low slope,
the acceleration phase (high slope), and dynamic rupture propagation with its speed comparable to
the shear wavespeed of the polycarbonate (Fig. S22d,e). The nucleation process is characterized
by the propagation of a pulse-like rupture front.

Figure S23 shows the evolution of rupture front position and rupture length in these sim-
ulations. The growth of the rupture length during the quasi-static phase can be fit by the same
functional form that we used in fitting the seismic observations, namely Eq. (3), plotted as ma-
genta dashed curves in Figure S23. The exponents *k* = 2.5–3.2 are similar to those estimated from
the high-frequency (HF) ground velocity envelopes in our observations (e.g., Fig. S16). The expo-
nent increases with effective normal stress because the rupture front is more quickly accelerating
for a larger effective normal stress. The flat line (at rupture length ~0.1 mm) in Figure S23 is an
artifact arising from the limit of spatial resolution; the rupture length is essentially zero over the
flat part. The duration of the quasi-static phase is shorter for smaller normal stresses because the
slip during the dynamic rupture becomes smaller, which leads to a shorter inter-event period and
duration of the quasi-static phase.

Guided by the results of numerical modeling, we interpret the HF seismic observations as
follows. The similarity of growth exponents of HF ground velocity envelopes (v) and slip zone
radius (R) is consistent with the following two different models. In the ‘ring model,’ the sources
of HF waves are uniformly distributed over a ring of finite width *W* near a circular slip front (HF area *A* = π*R*W) and their HF waves interfere coherently (v ~ A). In the 'disc model,' the HF sources are uniformly distributed over the whole slip area (A = π*R*²) and their HF waves interfere incoherently (v² ~ A). Both models lead to v ~ R and hence same growth exponent of velocity envelopes and slip zone radius. In both models the radiated HF power per unit of fault area is uniform, as if HF waves were generated by small foreshocks on a uniformly distributed set of small asperities.

Our numerical modeling results imply that the observed growth of high-frequency seismic
signals—either via the ring or disc models—is the nucleation phase of an earthquake rupture that
is characterized by the propagation of a slow-slip pulse. Very-low-frequency waves emerge from
both the HF sources and the slow slip.

**Data availability.** All seismic data used in this study, notably from the FLATS (XV) seismic network (doi:10.7914/SN/XV_2014), are publicly available from the IRIS Data Management Center (http://ds.iris.edu/ds/nodes/dmc/). Expanded results of moment tensor inversions are available in ref. ⁵⁶.

1. Peng, Z. & Gomberg, J. an integrated perspective of the continuum between earthquakes and slow-slip phenomena. *Nature Geoscience* **3**, 599–607 (2010).

2. Dodge, D. A., Beroza, G. C. & Ellsworth, W. L. Detailed observations of California foreshock sequences: Implications for the earthquake initiation process. *J. Geophys. Res.* **101**, 22,371–22,392 (1996).
3. McGuire, J. J., Boettcher, M. S. & Jordan, T. H. Foreshock sequences and short-term earthquake predictability on East Pacific Rise transform faults. *Nature* **434**, 457–461 (2005).

4. Bouchon, M. et al. Extended nucleation of the 1999 M<sub>w</sub> 7.6 Izmit earthquake. *Science* **331**, 877–880 (2011).

5. Bouchon, M., Durand, V., Marsan, D., Karabulut, H. & Schmittbuhl, J. The long precursory phase of most large interplate earthquakes. *Nature Geoscience* **6**, 299–302 (2013).

6. Kato, A., Fukuda, J., Kumazawa, T. & Nakagawa, S. Accelerated nucleation of the 2014 Iquique, Chile Mw 8.2 earthquake. *Sci. Rep.* **6**, 1–9 (2016).

7. Ogata, Y. Seismicity and geodetic anomalies in a wide area preceding the Niigata-Ken-Chuetsu earthquake of 23 October 2004, central Japan. *J. Geophys. Res.* **112** (2007).

8. Scuderi, M. M., Marone, C., Tinti, E., Di Stefano, G. & Collettini, C. Precursory changes in seismic velocity for the spectrum of earthquake failure modes. *Nature Geoscience* **9**, 695–702 (2016).

9. Schurr, B. et al. Gradual unlocking of plate boundary controlled initiation of the 2014 Iquique earthquake. *Nature* **512**, 299–302 (2014).

10. Mavrommatis, A. P., Segall, P., Uchida, N. & Johnson, K. M. Long-term acceleration of aseismic slip preceding the M<sub>w</sub> 9 Tohoku-oki earthquake: Constraints from repeating earthquakes. *Geophys. Res. Lett.* **42**, 9717–9725 (2015).
11. Obara, K. & Kato, A. Connecting slow earthquakes to huge earthquakes. *Science* **353**, 253–257 (2016).

12. Ide, S., Beroza, G. C., Shelly, D. R. & Uchide, T. A scaling law for slow earthquakes. *Nature* **447**, 76–79 (2007).

13. Kato, A. et al. Propagation of slow slip leading up to the 2011 M$_w$ 9.0 Tohoku-Oki earthquake. *Science* **335**, 705–708 (2012).

14. Campillo, M. & Ionescu, I. R. Initiation of antiplane shear instability under slip dependent friction. *J. Geophys. Res.* **102**, 20,363–20,371 (1997).

15. Dieterich, J. H. Earthquake nucleation on faults with rate- and state-dependent strength. *Tectonophysics* **211**, 115–134 (1992).

16. Lapusta, N. & Rice, J. R. Nucleation and early seismic propagation of small and large events in a crustal earthquake model. *J. Geophys. Res.* **108** (2003).

17. Rubin, A. M. & Ampuero, J.-P. Earthquake nucleation on (aging) rate and state faults. *J. Geophys. Res.* **110** (2005).

18. Kaneko, Y. & Ampuero, J.-P. A mechanism for preseismic steady rupture fronts observed in laboratory experiments. *Geophys. Res. Lett.* **38** (2011).

19. Ohnaka, M. & Shen, L. Scaling of the shear rupture process from nucleation to dynamic propagation: Implications of geometric irregularity of the rupturing surfaces. *J. Geophys. Res.* **104**, 817–844 (1999).
20. Ben-David, O., Cohen, G. & Fineberg, J. The dynamics of the onset of frictional slip. *Science* **330**, 211–214 (2010).

21. Nielsen, S., Taddeucci, J. & Vinciguerra, S. Experimental observation of stick-slip instability fronts. *Geophys. J. Int.* **180**, 697–702 (2010).

22. McLaskey, G. C. & Kilgore, B. D. Foreshocks during the nucleation of stick-slip instability. *J. Geophys. Res. Solid Earth* **118**, 2982–2997 (2013).

23. Latour, S., Schubnel, A., Nielsen, S., Madariaga, R. & Vinciguerra, S. Characterization of nucleation during laboratory earthquakes. *Geophys. Res. Lett.* **40**, 5064–5069 (2013).

24. Ellsworth, W. L. & Beroza, G. C. Seismic evidence for an earthquake nucleation phase. *Science* **268**, 851–855 (1995).

25. Mori, J. & Kanamori, H. Initial rupture of earthquakes in the 1995 Ridgecrest, California sequence. *Geophys. Res. Lett.* **23**, 2437–2440 (1996).

26. Ito, Y. et al. Episodic slow slip events in the Japan subduction zone before the 2011 Tohoku-Oki earthquake. *Tectonophysics* **600**, 14–26 (2013).

27. Linde, A. T., Suyehiro, K., Miura, S., Sacks, I. S. & Takagi, A. Episodic aseismic earthquake precursors. *Nature* **334**, 513–515 (1988).

28. Tape, C., West, M., Silwal, V. & Ruppert, N. Earthquake nucleation and triggering on an optimally oriented fault. *Earth Planet. Sci. Lett.* **363**, 231–241 (2013).
29. Linde, A. T., Gladwin, M. T., Johnston, M. J. S., Gwyther, R. S. & Bilham, R. G. A slow earthquake sequence on the San Andreas fault. *Nature* **383**, 65–68 (1996).

30. Ratchkovski, N. A. & Hansen, R. A. New constraints on tectonics of interior Alaska: Earthquake locations, source mechanisms, and stress regime. *Bull. Seismol. Soc. Am.* **92**, 998–1014 (2002).

31. Tape, C. et al. Transtensional tectonics of the Minto Flats fault zone and Nenana basin, central Alaska. *Bull. Seismol. Soc. Am.* **105**, 2081–2100 (2015).

32. Pollitz, F. F., Stein, R. S., Sevilgen, V. & Bürgmann, R. The 11 April 2012 east Indian Ocean earthquake triggered large aftershocks worldwide. *Nature* **490**, 250–253 (2012).

33. Ito, Y., Obara, K., Shiomi, K., Sekine, S. & Hirose, H. Slow earthquakes coincident with episodic tremors and slow slip events. *Science* **315**, 503–506 (2007).

34. Ghosh, A., Huesca-Pérez, E., Brodsky, E. & Ito, Y. Very low frequency earthquakes in Cascadia migrate with tremor. *Geophys. Res. Lett.* **42**, 3228–3232 (2015).

35. Ide, S., Imanishi, K., Yoshida, Y., Beroza, G. C. & Shelly, D. R. Bridging the gap between seismically and geodetically detected slow earthquakes. *Geophys. Res. Lett.* **35** (2008).

36. Ide, S. Characteristics of slow earthquakes in the very low frequency band: Application to the Cascadia subduction zone. *J. Geophys. Res. Solid Earth* **121**, 5942–5952 (2016).

37. Ide, S. A Brownian walk model for slow earthquakes. *Geophys. Res. Lett.* **35** (2008).
38. Gomberg, J., Agnew, D. C. & Schwartz, S. Y. Alternative source models of very low frequency events. *J. Geophys. Res. Solid Earth* **121**, 6722–6740 (2016).

39. Hawthorne, J. C. & Ampuero, J.-P. A phase coherence approach to identifying co-located earthquakes and tremor. *Geophys. J. Int.* **209**, 623–642 (2017).

40. Shearer, P. M., Prieto, G. A. & Hauksson, E. Comprehensive analysis of earthquake source spectra in southern California. *J. Geophys. Res.* **111** (2006).

41. Rogers, G. & Dragert, H. Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip. *Science* **300**, 1942–1943 (2003).

42. Bosl, W. J. & Nur, A. Aftershocks and pore fluid diffusion following the 1992 Landers earthquake. *J. Geophys. Res.* **107** (2002).

43. Renard, F. et al. Critical evolution of damage toward system-size failure in crystalline Rock. *J. Geophys. Res. Solid Earth* **123**, 1969–1986 (2018).

44. Kaneko, Y., Nielsen, S. B. & Carpenter, B. M. The onset of laboratory earthquakes explained by nucleating rupture on a rate-and-state fault. *J. Geophys. Res. Solid Earth* **121**, 6071–6091 (2016).

45. McLaskey, G. C. & Lockner, D. A. Preslip and cascade processes initiating laboratory stick slip. *J. Geophys. Res. Solid Earth* **119**, 6323–6336 (2014).

46. Dublanchet, P. The dynamics of earthquake precursors controlled by effective friction. *Geo-phys. J. Int.* **212**, 853–871 (2018).
47. Noda, H., Nakatani, M. & Hori, T. Large nucleation before large earthquakes is sometimes skipped due to cascade-up—Implications from a rate and state simulation of faults with hierarchical asperities. *J. Geophys. Res. Solid Earth* **118**, 2924–2952 (2013).

48. Rice, J. R. Spatio-temporal complexity of slip on a fault. *J. Geophys. Res.* **98**, 9885–9907 (1993).

49. Ye, L., Lay, T., Kanamori, H. & Rivera, L. Rupture characteristics of major and great ($M_w \ge 7.0$) megathrust earthquakes from 1990 to 2015: 2. Depth dependence. *J. Geophys. Res. Solid Earth* **121**, 845–863 (2016).

50. Zhu, L. & Helmberger, D. Advancement in source estimation techniques using broadband regional seismograms. *Bull. Seismol. Soc. Am.* **86**, 1634–1641 (1996).

51. Vallée, M. Stabilizing the empirical Green function analysis: Development of the projected Landweber method. *Bull. Seismol. Soc. Am.* **94**, 394–409 (2004).

52. Ni, S., Helmberger, D. & Kanamori, H. Energy radiation from the Sumatra earthquake. *Nature* **434**, 582 (2005).

53. Lapusta, N. & Liu, Y. Three-dimensional boundary integral modeling of spontaneous earthquake sequences and aseismic slip. *J. Geophys. Res.* **114** (2009).

54. Dieterich, J. H. Modeling of rock friction: 1. Experimental results and constitutive equations. *J. Geophys. Res.* **84**, 2161–2168 (1979).
55. Ruina, A. L. Slip instability and state variable friction laws. *J. Geophys. Res.* **88**, 10,359–10,370 (1983).

56. Silwal, V. Seismic moment tensors for six events in the Minto Flats fault zone, 2012–2016 (2018). ScholarWorks@UA at http://hdl.handle.net/11122/8253: descriptor file, text file of catalog, figures with waveform fits, and input weight files.

**Corresponding author** Correspondence and requests for materials should be addressed to Carl Tape.

**Acknowledgements** Seismic instruments, data archiving, and data access were supported by NSF grants EAR-1352668 and EAR-1249780, the Alaska Earthquake Center, the IRIS Data Management Center, and the PASSCAL Instrument Center.

**Author contributions** All authors contributed to the manuscript. Tape, Holtkamp, and Silwal performed the seismological analyses of the two events. Hawthorne performed the phase coherence analysis (Fig. 3). Kaneko, Ampuero, and Hawthorne contributed the rate-and-state modeling and interpretation. Ji performed the source time function estimation. Tape, Ruppert, Smith, and West were responsible for the deployment of FLATS seismic stations. Ruppert and Holtkamp discovered the high-frequency and low-frequency signals for the 2015 event. Holtkamp discovered the low-frequency foreshock (VLFE) of the 2016 event.

**Additional information** Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints.

**Competing financial interests** The authors declare that they have no competing financial interests.
<figure><img src="image_2.png" /><figcaption>Figure 1: Minto Flats fault zone of central Alaska. (a) Beachballs represent source mechanisms for earthquakes (EQ) and very-low-frequency earthquakes (VLFE) discussed in this study (Tab. S1). From north to south, these are: 1995-10-06 $M_w$ 6.0 EQ (18 km), 2015-09-12 $M_w$ 3.8 VLFE (21 km), 2012-04-11 $M_w$ 3.8 VLFE+EQ (16 km), 2015-10-22 $M_w$ 2.6 EQ (18 km), 2013-03-12 $M_w$ 3.5 VLFE (23 km), 2016-01-14 $M_w$ 3.7 VLFE+EQ (17 km), 2015-10-31 $M_w$ 3.4 EQ (25 km).</figcaption></figure>

The thick black lines denote the two left-lateral faults of the Minto Flats fault zone. Inverted triangles denote seismic stations, with MDM at far right. (b) Envelopes of high-frequency seismograms at station MDM for five events in (a). For a typical earthquake (top), there is no high-frequency foreshock (HFF) signal prior to the P wave of the earthquake. Other events show an increasing signal that is associated with nucleation and VLFEs. The arrows labeled VLFE indicate the onset time of low-frequency waves such as those shown in Figure 2a.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Figure 2: Seismic waveforms for the 2015 and 2016 events observed at stations FTGH (red traces) and F6TP (black traces). (a) The 2015 very-low-frequency earthquake (VLFE). Top: seismogram causal-filtered 20–100 s. Middle: seismogram filtered 2–8 Hz. Bottom: envelope of the high-frequency seismogram, plotted with log scaling. (b) The 2016 VLFE+Earthquake, with high-frequency foreshock (HFF) nucleation signals. Top: seismogram causal-filtered 20–100 s, cut at the earthquake S time and zoomed in on the vertical scale to show the shape of the waveform, with an expanded view of the precursory VLFE in the gray box above. Middle: seismogram filtered 2–8 Hz. Bottom: envelope of the high-frequency seismogram. The dashed line is the P wave arrival for the *M*<sub>w</sub> 3.7 earthquake.
<figure><img src="image_5.png" /><figcaption>Figure 3: Phase coherence between foreshock signal and mainshock signal for the 2016 event, filtered between 1–10 Hz. Each subplot is a 4 s time window in a north-south oriented vertical plane; time 0 s is the mainshock P wave.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Alternative interpretations for the two stages of nucleation and rupture. Our preferred interpretation for the 2016 event is **a→b**, that the *M*<sub>w</sub> 3.7 earthquake nucleated from a very-low-frequency earthquake (VLFE) associated with slow slip. Alternative interpretations are **a→d**, **c→b**, and **c→d**. The 2015 VLFE event comprised stage 1 (a or c) but did not transition into an earthquake.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

Figure 5: Comparison between rate-state modeling and observations of power-law growth of the rupture stage interpreted as nucleation. (a) Evolution of non-dimensional rupture length as a function of non-dimensional time for the cases shown in Figure S23. The rupture length and time are non-dimensionalized by $h_{RR}^*$ and $\sigma(b-a)/\dot{\tau}$, respectively, where $h_{RR}^*$ is theoretical estimate of nucleation size, $\sigma$ is the effective normal stress, $b-a$ is the rate-and-state parameter, and $\dot{\tau}$ is the background loading rate$^{44}$. Time $t=0$ is defined as the time at the end of the acceleration phase.

(b) Dimensional rupture length vs. time. We assume that $h_{RR}^* = 100$ m and $\sigma(b-a)/\dot{\tau} = 0.15$ s.

(c) Example of the observed high-frequency foreshock (HFF) amplitude growth for the 2016 event for station MDM (see Figure 1b). See Figure S16 for more examples.