Rapid Tremor Migration During Few Minute-Long Slow Earthquakes in Cascadia

few minutes and a few hours, have not yet been observed geodetically but are frequently suggested by varying tremor migration and amplitude. In Cascadia, tremor often migrates 40–60 km along dip during hour-long tremor streaks, with migration rates 50–500 times faster than the Abstract Slow earthquakes are now commonly found to display a wide range of durations, moments, and slip and propagation speeds. But not all types of slow earthquakes have been examined in detail. Here we probe tremor bursts with durations between 1 and 30 min, which are likely driven by few minute-long bursts of aseismic slip. We use a coherence-based technique to detect thousands of tremor bursts beneath Vancouver Island in Cascadia. Then we examine 17 of the ruptures by tracking their evolving tremor locations over an 8-km region. We find that tremor migrates at rates of 3–25 m/s: faster than longer tremor bursts. Though some observational biases persist, the short events' speeds appear to fill a gap in the spectrum of observed slow earthquakes. They may provide further evidence that whatever fault zone process creates slow earthquakes, it must allow for faster slip and propagation in smaller ruptures.

Tremor migration has not yet been observed in detail on timescales shorter than 10 min, but several features of tremor suggest that complex, rapid propagation should continue to short timescales. First, tremor varies in amplitude on a range of timescales, from seconds to days (Ghosh et al., 2010;Ide, 2010a;Obara, 2002;Shelly et al., 2006), and those variations are correlated with aseismic deformation (W. Frank, 2016;Hawthorne & Bartlow, 2018;Hawthorne & Rubin, 2013b). Second, some 20-200 s-long increases in tremor amplitude are associated with 20-200 s-long increases in slow slip moment rate. These moment rate increases are observable in long-period seismic data, and the 20-50 s-long events are called very low frequency earthquakes (VLFEs) (e.g., Baba et al., 2020;Hutchison & Ghosh, 2016;Ito & Obara, 2006;Maury et al., 2016;Takeo et al., 2010;Walter et al., 2013).
In this study, we identify and analyze tremor bursts with durations between 1 and 30 min. Many of these events are slightly longer than VLFEs but shorter than previously detected (>10-min) tremor migrations; few-minute tremor bursts have not been analyzed in detail in previous work. To fill this observational gap in the slow earthquake spectrum, we first identify thousands of tremor bursts in Cascadia and then probe 17 of them in more detail. We identify rapid tremor migration which likely reflects rapid migration of aseismic slip.

Motivation to Constrain the Spectrum of Subevents
We analyze migration on few-minute timescales for two reasons: (a) because we wish to more fully observe the range of behaviors in slow earthquakes and (b) because the range of slip speeds and behaviors could help us determine which fault zone process creates slow earthquakes. The propagation speeds of few hour-long subevents have already been used to test some models of slow earthquakes (Ariyoshi et al., 2009;Luo & Ampuero, 2017;Luo & Liu, 2021;Rubin, 2011). Some researchers have modeled the range of slow earthquakes as parts of a diffusive process, where ruptures grow more quickly when they are small (Ando et al., 2012;Ide, 2008;Ide & Maury, 2018). Other researchers have proposed that rapid subevents could reflect the rupture of asperities or asperity clusters embedded in the slow slip region. They have successfully produced the propagation speeds of few hour-subevents (RTRs) with a relatively simple approach: by mixing unstable patches into a region with a nominally stable slow slip rheology. The unstable patches effectively increase the local stress drop and drive more rapid slip (Ando et al., 2012;Ariyoshi et al., 2009;Colella et al., 2012;Luo & Liu, 2021;Nakata et al., 2011;Peng & Rubin, 2018).
However, it may not be plausible that slip speed increases by many orders of magnitude simply because the stress drop that drives rupture increases. For many of the proposed slow slip rheologies, the stress drop required for rupture increases dramatically as the rupture speed increases (e.g., Hawthorne & Rubin, 2013c). It may be that increased local stress drops can provide enough energy to increase the rupture speed by a factor of 10-20, as seen in RTRs, but a factor of 100 or 1,000 rupture speed increases may require a spatial variation in the resistance to accelerating slip (e.g., Hawthorne et al., 2016).
Such spatially variable resistance to slip is intriguing because it should not exist for several of the proposed slow slip rheologies. For instance, if slow slip events happen because the rheology at depth imposes a temperature-dependent speed limit (e.g., Hawthorne & Rubin, 2013a;Shibazaki & Iio, 2003;Shibazaki & Shimamoto, 2007), that speed limit should stay roughly the same throughout the slow slip region, where the temperature stays relatively uniform.
There is, of course, already evidence that slow earthquake slip rates vary by at least four orders of magnitude, from 10 −7 m/s in slow slip to 0.1 or 1 mm/s in VLFEs and tremor (e.g., Bartlow et al., 2011;Bostock et al., 2015;Dragert et al., 2001). However, it remains controversial whether tremor, VLFEs, and tremor bursts are created by the same rheology that governs slow slip rates. A single fault zone rheology is suggested by a systematic trend in observed slow earthquakes: smaller events are faster. In Figure 1b, we plot the propagation speeds and durations of tremor bursts from a variety of studies, and we see that shorter tremor bursts migrate faster. In Figure 1c, we 3 of 19 plot the moments and durations of slow earthquakes from a variety of studies, as was done by Ide et al. (2007). The data used in the plot are listed in Table S3. From these data, we see that observed slow earthquakes' moments M 0 scale roughly linearly with their duration T.
It will be useful to understand the slip and propagation rates implied by such a linear moment-duration scaling. To first order, moment in a slip event is proportional to slip times area: to δR 2 , where δ is the spatially averaged slip and R is the radius or long axis of the slip event. But slip δ is proportional to ΔτR: to the stress drop times the rupture radius. If we assume that slow earthquakes have magnitude-independent stress Horizontal bars indicate uncertainties in propagation velocity given a 0.5-km uncertainty in tremor location. The dashed parallelogram bounds the types of events we can observe with our approach. The solid black lines in panels (a-c) indicate a propagation rate that scales as T −2/3 and a moment equal to 3 × 10 12 N m s −1 times T. We map the line from panel c to the panels a and b assuming an elliptical rupture with a 3:1 aspect ratio and a 30-kPa stress drop. (b) Duration versus propagation velocity and (c) duration versus moment for our observations as well as for a selection of previous studies, indexed by the numbers below. Note that trends are visible in some studies but that there is often more uncertainty when comparing between locations. To avoid clutter, we plot only a handful of observations randomly selected from each study when a large number of events are detected. Many authors publish only a single average propagation rate and uncertainty, or they plot a handful of figures. In those cases, we choose one or a few values from the published distribution or extract rough propagation rates from the figures. Values are taken from 1: Shelly (2017) drops, as weakly suggested by observations of slow slip events, RTRs, and low frequency earthquakes (LFEs) (Chestler & Creager, 2017;Gao et al., 2012;Hawthorne et al., 2016), we can rewrite moment in several ways. First, 0 ∝ 2 ∝ Δ 3 ∝ ( ) 3 , and second, 0 ∝ 2 ∝ Δ −2 3 ∝ (̇) 3 . Here we have inserted an event-average propagation rate V p ∝ R/T and an event-averaged slip rate ̇∝ ∕ . If we insert M 0 ∝ T into these equations and solve for V p and ̇ , we find that a linear moment-duration scaling coupled with a moment-independent stress drop implies that both propagation rate V p and slip rate ̇ scale as −2∕3 0 . Since the linear moment-duration trend appears to extend all the way from M w 7 slow slip events to M w 1 LFEs, over a factor of 10 9 change in moment, we may expect a factor of 10 6 change in slip rate. It would seem sensible to start assessing which rheologies would allow slip speeds that are 10 4 times faster on 400-m LFE patches than on 400-km slow slip regions.
At this point, however, it is also sensible to recall that the scalings between moment, duration, and slip rate remain uncertain. Other moment-duration scalings have been observed. Bostock et al. (2015) and Farge et al. (2020) inferred a very weak moment-duration scaling, with ∼ 0 0 or 0.1 0 , among the LFEs that compose tremor. Michel et al. (2019) and Supino et al. (2020) identified a moment-duration scaling similar to that seen in normal earthquakes, with ∼ 1∕3 0 among geodetically observed slow slip events and among sub-second tremor bursts, respectively. Further, the linear moment-duration trend identified by Ide et al. (2007) depends on connecting days to months-long slow slip events and seconds-long LFEs and VLFEs ( Figure 1c). There are significant gaps along that trend. We now have abundant observations of sub-second LFEs, 10-s VLFEs, and days to months-long slow slip events, but there are fewer observations between those durations, and there may be missing observations that fall off the trend (Gomberg et al., 2016).
To overcome all the observational gaps, we will require a wide range of approaches, with a number of studies. For instance, recent work has suggested an expansion of the tremor band to longer durations Masuda et al., 2020). But here we focus on a different band and attempt to expand the range of slow slip subevents to shorter durations: between 2 and 10 min. We modify our tremor detection methodology to search for the expected short, rapid migration on these timescales, and we partly fill the apparent gap in the duration-propagation rate trend (Figure 1b). We note, however, that our approach still suffers from observational bias; it was not designed to find events that are much faster or slower than the along-trend speeds.
In the sections that follow, we first describe our phase coherence-based tremor detection approach (Section 3) and the available data and processing (Section 4). Then we describe the observed large-scale tremor migration patterns in Section 5, examine small-scale tremor migration in 17 tremor bursts in Section 6, and discuss the migrations' implications in Section 7.

Tremor Detection Method: Identifying Tremor Near Template Locations
To identify and locate tremor, we use a phase coherence-based approach developed by Hawthorne and Ampuero (2017), which is a variant of empirical matched field techniques (e.g., Bucker, 1976;Corciulo et al., 2012;Harris & Kvaerna, 2010;Wang et al., 2015). This approach allows us to identify tremor that ruptures fault patches close to known LFEs. The coherence calculation is able to identify tremor even if the tremor ruptures are complex or if the tremor consists of a series of ruptures, as the method combines two common approaches to identifying tremor. First, as inspired by matched filter techniques (Bostock et al., 2012;Brown et al., 2008;W. B. Frank et al., 2014;Shelly, 2017), the calculation compares seismograms between events. It assesses whether the template and target signals could have the same Green's functions: if they result from the same source-station path. Second, as inspired by cross-station techniques (Armbruster et al., 2014;Peng et al., 2015;, the calculation compares seismograms between stations or components. It assesses whether the signals at all stations or components could result from the same tremor source time functions.

Inter-Station Coherence
In all of the coherence calculations, we begin with a set of template seismograms d tkm that were created by Bostock et al. (2012) and which represent the signals generated by LFEs occurring at 130 locations on the plate interface (red crosses in Figure 2). The LFEs are recorded at a range of stations k and on three components m (east, north, and up). We compare the template seismograms at each station with 30-60-s-long intervals of target seismic data (d dkm ). To assess whether a 30 or 60-s interval contains tremor coming from the same location as the template, we compute the inter-station phase coherence at a range of frequencies f: Here ̂( ) and ̂( ) are the Fourier transforms of the template and target data at station k and component m, and we compare between stations k and l in each term. There are N stations in total, and we average over the N(N − 1)/2 station pairs and over the three components at each station. Note that we have dropped the frequency indexing on the right hand side for readability, and we also average over frequencies f between 1 and 6 Hz.
Note that if the target seismograms d dkm record tremor from the same location as the template seismograms, then the template and target seismograms may be written as ̂=̂̂ and ̂=̂̂ , where and are the template and target tremor source time functions, is the path effect, and becomes So by identifying intervals with high phase coherence , near 1, we can identify intervals when tremor is occurring at the same location as previously located templates. Synthetic tests suggest that is significantly larger than zero only when tremor occurs within about 0.5 km from the template: within a fraction of one seismic wavelength ( Figure S3 in Supporting Information S1). We will therefore use this location-specific calculation in Section 6 as we track the spatial evolution of tremor locations.

Inter-Component Coherence
In some cases, however, we do not need or want 0.5-km precision. For instance, in Section 5, we will identify bursts of tremor. In that case, we want to identify tremor that is roughly the same area as the template: within 10 km or so. In such situations, we compute an inter-component phase coherence: Here we multiply the Fourier domain seismograms across components m and n rather than across stations k. Then we average over component pairs and over stations.
is often high (closer to 1) even when the target tremor is slightly offset from the template. To understand why, note that when the tremor is close to the template, its Green's functions are likely to have shapes similar to the template's Green's function. The tremor Green's functions ′ ( ) are approximately time-shifted versions of the template Green's functions g km (t). They may be approximated by ′ = ( − Δ ) , where the time shift Δt results from the difference in travel time from the source. Now we may note that if we have multiple recordings on the same station, just at different components m and n, the change in travel time Δt will remain the same. If we input tremor with these shifted Green's functions into the phase coherence calculation in Equation 3, we eliminate the travel time change and obtain Here s t and s d are the source time functions of the template and tremor signals, respectively.
It is of course difficult to know whether the Green's functions' shapes remain the same over broad regions. We find empirically that the Green's functions retain a similar enough shape for detection even as tremor locations change by 10-20 km; we sometimes obtain high for one template when the inter-station calculation is high only for templates 10 or 20 km away.
To assess spatial coherence somewhat more systematically, in Figure 3 and Figure S4 in Supporting Information S1 we compare the times of high coherence at templates with varying distances. We first identify times when coherence is more than 2σ ( Figure S4 in Supporting Information S1) or 2.5σ ( Figure 3) above the background at one template. Then we determine the fraction of those times when coherence is more than 2 or 2.5σ above background at a neighboring template. We subtract the fraction of the time expected for chance detections. The consistency is never perfect; even templates 1 km apart identify simultaneous high coherence (with a 2.5σ threshold) only 30% of the time. However, the simultaneous detection rate decays slowly, reaching 10% when the templates are more than 10 km apart. Given that tremor may be generated on either side of the templates, not just between them, the slow decay suggests that inter-component coherence can detect tremor over 10-20-km-wide regions.

Templates, Data, and Processing
In our initial approach to the data, we compute the phase coherence and between each template and the seismic data recorded during 13-20 days-long intervals during four major slow slip events in 2004, 2008, 2009, and 2010. The 130 LFE templates created by Bostock et al. (2012) are located beneath the southern tip of Vancouver Island and the Juan de Fuca Strait, at depths ranging from 28 to 45 km (crosses in Figure 2).
We use data from stations in permanent and temporary seismic networks, including the Canadian National Seismograph Network (CNSN), the POLARIS-BC Network (PO; Nicholson et al., 2005), the Plate Boundary Observatory Borehole Seismic Network (PBO), the Pacific Northwest Seismic Network (PNSN), and the USArray Transportable Array (TA). The available networks and stations evolved between the different slow slip events. Stations from the POLARIS-BC network were available only during the 2004 event while the PBO stations are available only after 2008. The POLARIS-BC network is ideal for this study, thanks to its dense configuration across the southern half of Vancouver Island, and Figure 2 shows the stations used for the 2004 slow slip event. Maps and tables of the 2008, 2009, and 2010 networks are available in Figure S1 and Table S2 of Supporting Information S1.
To prepare the data, we filter the target and template seismograms to between 0.6 and 20 Hz and downsample to a common sampling rate of 40 Hz. We extract a portion of each template, from 0.2 s before to 4.8 s after a manually picked S-wave arrival and then compute the coherences and between these template segments and a set of overlapping 60 s-long windows of target data, starting every 6 s.
Further details of the C p calculations are as described by Hawthorne and Ampuero (2017), but to summarize, we first cross-correlate at each station in the time domain, computing d dkm ⋅d tkm . Then we taper the time domain correlation with a Hanning filter, convert to the frequency domain, compute C p , and average C p over frequencies between 1 and 6 Hz. When the main slow slip front reaches the locations of the templates, the coherence values begin to spike to values well above noise-induced variability. The templates in Figure 4 see the main front arrival on the 14th and 15th July, respectively, and have values that reach 0.38 and values that reach 0.17. However, the coherence values are not continuously high. The activity is fragmented into 1-30 min-long spikes separated by intervals of low coherence that last minutes to hours. The spikes associated with the templates in Figure 4 can be seen in more detail in the 8 hr-long windows in panels c and d, but similar spikes in coherence are observed for all 130 templates examined in this study. The most intense sequence of spikes typically lasts 1-2 days, while the main front passes, but spikes in coherence can be seen for up to 5 days.

Observed Large-Scale Tremor Patterns
The spikes in the phase coherence and are presumably created by bursts of tremor occurring on the plate interface. The inter-component coherence is ideal for identifying and measuring the duration of these bursts, as seems to remain high even when tremor spreads to locations as far as 10-20 km from a given template. In contrast, the inter-station coherence decreases when the tremor is offset by more than a fraction of the seismic wavelength. The difference in inter-component and inter-station tremor detection is apparent for a number of spikes for template #246 (Figures 4a and 4c). For instance, the spike at 3:45 on 15 July is longer on the time series, and the spike at 4:30 on 15 July is observable only the time series.
We use a simple peak detection algorithm included in SciPy to identify a number of spikes in each time series. Proposed spikes are identified as maxima in when where the standard deviation of the phase coherence during a few-day interval before tremor begins each year, and α is a factor between 2.0 and 3.5. The beginning and end of the spikes are the times when decreases to half its maximum value, and we accept spikes that last at least 60 s. Depending on the threshold α, we detect between 22,000 and 86,000 events in 2004. Several spikes are delineated in yellow in Figure 5a. Catalogs of the spikes are provided as a Supporting Information S1, and Figure S7 in Supporting Information S1 shows the distribution of spike durations for different values of α. Note that because we use a 60-s window for our coherence calculations, is smoothed on that timescale, and the durations of shorter spikes may be overestimated. Figure 2 shows the spatio temporal distribution of tremor bursts detected with α = 3.0. Each burst is plotted at a random location within 5 km of the associated template location and are colored by time. The bursts track the along-strike propagation of the main ETS front in 2004, from the Juan de Fuca Strait on 10 July to ∼90 km northwest of the Strain on 19 July (e.g., Bostock et al., 2015;Wech et al., 2009).
We can also compare our results directly with the LFE detections of Bostock et al. (2015). Vertical red lines in Figure 5 mark LFE detections with the relevant template. The matched filter detections in the catalog coincide remarkably well with times of high phase coherence. All LFE detections occur within intervals of high , though a few intervals of high do not include a LFE detection. The lack of LFE detections in some high intervals may arise because tremor is coming from an adjacent part of the fault or because the tremor time series is complex, so that it is difficult to separate overlapping LFEs with a matched filter approach.
We observe similar tremor burst spacing and migration patterns for the 2008, 2009, and 2010 slow slip events ( Figure S1 in Supporting Information S1). 8 hr-long time series from the four events can be compared in Figure 5, though we focus on the 2004 results in this study because the C p time series have the highest resolution in 2004, when the POLARIS seismic network was running on Vancouver Island.

One Example
Our observed spikes in , along with tremor spikes seen in previous work Ghosh et al., 2009;Rubin & Armbruster, 2013), suggest that much of the tremor in Cascadia occurs in short bursts. Here we seek to probe the bursts in more detail. We examine the shape and migration of some of the shorter bursts.
We track the spatial and temporal evolution of 17 tremor bursts that are visible as well-resolved spikes in the record. We identify a high-quality template that records each burst and then define a circular grid of potential tremor locations around that template, as illustrated in Figure 6a and Figures S8-S24 in Supporting Information S1. Each grid is 8 km in diameter and is inclined along the slab interface identified by McCrory et al. (2012). In order to track tremor within the grid, we note that tremor coming from each of the possible locations is likely to have a Green's function whose shape is similar to the template's Green's function. We verify that similarity by comparing the waveforms of closely spaced template LFEs; the waveforms of templates located about 5 km apart have similar shapes (see Figure S2 in Supporting Information S1).
We shift the timing of the template seismograms to reflect the variation in the source-station travel time among the grid locations. The travel time for each location is computed using a uniform shear wave velocity model. We have estimated the apparent shear wave velocity for each LFE template by plotting the variation in 3-D distance from the LFE to the various stations against the arrival time for each station. We observe a linear relationship between distance and travel time, suggesting that a uniform velocity model is sufficient for our analysis. Tests with a layered velocity model and ray path calculations achieved similar results, presumably because the relative location shifts give similar time shifts for these velocity models.
Once we have time-shifted the template waveforms for a given location, we compute the inter-station phase coherence to determine when tremor occurs at that location. We compute for each point on the grid, using 1 min windows spaced every 6 s, so with a large overlap. Then we visually examine the patterns in to identify any migration. As one example, Figure 6 shows snapshots of the coherence during one three minute-long burst. During this time, the region of high phase coherence migrates about 1.6 km at a speed around 33 km/hr. The tremor moves from southeast to northwest, roughly along the strike of the subduction zone. This northwestward migration is pulse-like; the first location stops generating tremor before the last location generates tremor.

Propagation for 17 Tremor Bursts
Tremor migration is also well-resolved for 16 other analyzed bursts, with migration durations between 60 and 1,100 s. These are illustrated in Figures S8-S24 and flipbooks M1-M6 in Supporting Information S1. To more precisely characterize the tremor migration speed in each burst, we visually estimate the propagation azimuth and then create profiles of coherence along that azimuth and along other azimuths within 10 or 20°, as seen during one rupture in Figure 7. For the 17 bursts examined, we find that the rupture front can be reasonably approximated by the location where first exceeds 0.02, so we compute a linear regression between this front position and time. The slope gives the propagation velocity along the proposed azimuths ( Figure S6 in Supporting Information S1). We take the preferred propagation azimuth to be that with the fastest propagation velocity.
The 17 analyzed bursts and their propagation speeds are listed in Table S2 of Supporting Information S1. Figure 1a shows the relationship between the bursts' duration T and the propagation velocity V r of the 17 events. Horizontal bars in Figure 1a give rough uncertainties in the propagation velocities. To obtain these bounds, computed for each point on the grid shown in a and interpolated. Indicated times are the middle of one minute-long windows used to compute the phase coherence. Red stars mark the location of the low frequency earthquake used as template. In this example, the slow-slip propagates roughly 1.6 km at 9.2 m/s, or 33 km/hr. we allow 0.5-km uncertainties in the starting and ending locations of the propagating fronts, as implied by the synthetic calculations in the next section.
The observed duration-propagation velocity relationship may or may not be representative of the general population of tremor bursts. We did not choose the bursts to be analyzed in a systematic way, as the main goal of this study is simply to look for rupture propagation in few minute tremor bursts. We simply selected the 100 tremor bursts with the highest maximum values, tracked the spatial evolution of in each of the 100 bursts, visually identified the bursts with clear migration, and probed only those events in more detail. The remaining tremor bursts, which do not show clear migration, may move too slowly or too quickly for us to see, or the tremor may come from outside the 8-km circle where we compute .
While we should keep in mind that we do not select our analyzed bursts very rigorously, it remains interesting to note that among the 17 events analyzed, shorter events propagate faster. The minute-long events propagate at more than 20 m/s while the 15 min-long events propagate at only ∼4 m/s. The faster propagation of shorter bursts persists for two definitions of burst duration. The open squares in Figure 1a indicate the durations of visually identified tremor migration while the filled circles indicate durations estimated from the time series: when is above a local background value. This latter definition of duration is likely to be more accurate, as the tremor could migrate out of the 8-km grid where we computed , and the inter-component phase coherence can identify at least some tremor in a broader region.

Range of Observable Propagation Velocities
However, before we interpret the observed durations and propagation velocities, it is important to note that we have chosen our methodology to examine a particular range of tremor bursts: those with durations between 1 and 30 min and propagation extents between 1 and about 8 km. This range of observable speeds and durations is outlined with a dashed line in Figure 1a.
We cannot identify migration over distances less than 1 km because of the resolution of the phase coherence calculations given 1-6 Hz seismic data from 5 to 10 stations. To map the potential smearing of the high coherence region for a given template, we follow an approach common in array analysis (Hawthorne & Ampuero, 2017;Rost & Thomas, 2002). We create a synthetic signal for each point on the circular grid by time-shifting the template waveforms. Then we compute between those synthetics and the unshifted template seismograms. Some of the noisier template seismograms allow for elongate smearing of the high over 3 km-long distances (Figures S3c and S3d in Supporting Information S1). However, we have looked for tremor propagation only with the higher-quality templates, which show relatively circular smearing of high over smaller regions, with half-width of around 0.5 km (Figures S3a and S3b in Supporting Information S1). These resolution tests imply that we should be able to identify tremor that propagates around 1 km or more.  Table  S2 in Supporting Information S1). The profile crosses the center of the grid and is orientated along the propagation direction. The dashed line indicates the threshold used to estimate the rupture velocity.
We also use synthetic tests to check that we can resolve tremor locations even if tremor moves slowly or quickly. We use a template LFE to generate seismograms from two tremor sources, moving at 4 m/s and 40 m/s, respectively. Then we analyze the seismograms as we would an observed tremor burst, and we successfully resolve 0.5 km of motion over 2 min, with a best-fitting rate of 4 m/s, as well as 4 km of motion over 2 min, with a best-fitting rate of 39 m/s ( Figure S6 in Supporting Information S1).
But even if we can identify slow and rapid migration, we cannot identify migration over large distances. For instance, we cannot track migration over distances larger than 8 km because that migration would extend outside the 8 km-wide circle where we look for tremor. We have not developed the technique to map migration from the area around one LFE template to another. In fact, all of the migration we observe stays within 2-3 km of the center of the circle; it may be that the Green's functions change or the inter-station time shifts are inaccurate at larger distances, so perhaps we can map migration only over distances less than 4-6 km. This migration tracking limit constrains the open squares to fall well below the upper (8-km) diagonal dashed line in Figure 1a.
However, we do not have to track migration over the entirety of a tremor burst in order to determine a migration speed; we can track migration during just part of the rupture. And we have another approach to infer rupture durations: the duration of the inter-component coherence . Comparison of detections between templates suggest that can sometimes detect tremor 10-20 km away (see Section 3.2). The filled circles in Figure 1a can thus in principle plot up to a factor of 2 above the 8-km line, as they indicate durations taken from spikes in the record. However, the true detection capability and thus the maximum duration of likely depends on the signal to noise ratio during the tremor bursts.
The roughly 1-8 km constraints suggest that we would observe an anticorrelation between burst duration and propagation speed even for a collection of bursts with random properties. However, the durations and propagation speeds we observe do not seem to fill the box of observable values; they are more consistent with the T 2−/3 trend that one would expect for slow earthquakes whose moments scale linearly with duration. It thus seems likely that our observed duration-propagation speed anticorrelation is real-not entirely an observational artifact, but since we did not choose the 17 bursts to analyze rigorously, we cannot be sure.

Discussion
We have used a high-precision, coherence-based technique to identify numerous bursts of tremor. We mapped tremor migration over 1-6 km in 17 bursts with durations between 1 and 22 min. The tremor migrates at speeds of 3-25 m/s, moving more quickly in shorter bursts. The sub-10 min, rapidly migrating bursts represent a new observation. They are yet another category of slow earthquakes that must be reproduced by any complete physical model of subduction zone slip.

Pulse-Like Ruptures
If we assume, as seems plausible, that the observed migration of tremor results from a migrating location of aseismic slip, it is interesting to note that the propagation of tremor is pulse-like rather than crack-like; the locations slipping early in the bursts (e.g., the brighter portions of Figures 6b and 6c) stop slipping before slip occurs at later locations (e.g., in Figures 6g and 6h). Such pulse-like migration of tremor and slip is also apparent in the large slow slip events in Cascadia, as well as in some longer tremor bursts (Dragert et al., 2001;Ghosh et al., 2010;Royer et al., 2015;Wech et al., 2009).
Pulse-like ruptures can appear unintuitive because slip at later locations should increase the stress at the initial locations, and that stress increase has the potential to drive slip. The pulse-like ruptures could indicate that the slow slip region has a particular type of rheology: one that allows a rapid recovery in stress as the slip rate slows, so that the initial location can accommodate an increasing stress as it slows down but other parts of the fault accelerate (Bizzarri, 2010;Heaton, 1990;Lu et al., 2007;Noda et al., 2009;Zheng & Rice, 1998). Alternatively, the pulse-like ruptures could indicate that the subevent rupture is constrained to an elongate region. Slip may migrate along the long axis of that region, and the initial slipping location may stop slipping because it has slipped enough relative to its short-axis edges to accommodate the local stress drop. Slip at the far end of the rupture may produce an insignificant stress change at the initial location (e.g., Dal Zilio et al., 2020;Hawthorne & Rubin, 2013a;Michel et al., 2017).
However, pulse-like ruptures are unlikely in some models of slow slip. One of the first proposed explanations of slow slip suggests that slow slip regions have a "standard," potentially unstable, velocity-weakening rheology but that the regions have a particular size; they may be large enough to accelerate but too small to reach seismic slip speeds (Li et al., 2018;Y. J. Liu & Rice, 2005Romanet et al., 2018;Rubin, 2008). But most simulations of those size-limited slow slip ruptures appear more crack-like than pulse-like, at least in a visual inspection (Y. J. Liu & Rice, 2005;Rubin, 2008). A more rigorous investigation of the slip rate profiles in these models would help us further assess whether fault sizes and geometry alone can explain slow slip events.

Too Fast to Be Driven by a Change in Stress Drop?
We may also investigate the rheology of the slow slip region by addressing their speed: how can 5 min-long subevents propagate at speeds of 10 m/s: 200 times faster than the main slow slip front, which moves around 5 km/day (Dragert & Wang, 2011;Wech et al., 2009)? Do the subevents have a greater driving stress drop and thus a larger strain energy release, or do they have a lower resistance to acceleration: a smaller fracture energy? To partially address this question, we may note that the strain energy released in an elongate rupture normally scales as Δτ 2 W: as the stress drop Δτ squared times the rupture width W (e.g., Lawn, 1993). This strain energy must equal the fracture energy dissipated by the rupture, which is a function of the rheology, the initial conditions, and the slip rate. Most of the complex rheologies proposed to explain slow slip have fracture energies that increase dramatically as the slip rate increases. The strong increase in resistance with slip rate keeps the slip rates low.
Our events have propagation speeds around 200 times faster than the main slip front. In models of propagating ruptures, slip rate is proportional to rupture speed, to first order (e.g., Ida & Aki, 1972;Rubin, 2011), so we may infer a slip rate 200 times faster than the slip rate in the main front. A factor of 200 increase in slip rate is likely to require at least a factor of 10 increase in fracture energy in a shear-induced dilatancy model (L. Liu et al., 2010;Segall et al., 2010) and at least a factor of 5 increase in fracture energy in a model with a velocity-strengthening transition (Hawthorne & Rubin, 2013a, 2013c. Our subevents have widths W at least a factor of 10 narrower than the ∼60-km wide main slow slip region, so for their slip to supply a factor of >5 increase in energy (for Δτ 2 W to go up by a factor of 5 or more), they would need stress drops at least 7 times larger than the main event stress drop. We do not have stress drop estimates for our subevents, but geodetic and tremor count-based analyses for half-to few-hour events suggest that subevent stress drops are comparable to or smaller than the main event stress drop (Bletery et al., 2017;Hawthorne et al., 2016;Rubin & Armbruster, 2013).
Some modelers have produced locally high stress drops and slip rates by mixing unstable patches into a mostly stable slow slip region (Ariyoshi et al., 2009(Ariyoshi et al., , 2012Colella et al., 2012;Luo & Liu, 2021;Peng & Rubin, 2018). However, these patch-driven models have focused on slower tremor fronts than those examined here. They have so far allowed propagation rates less than 10-50 times faster than the main front (Ariyoshi et al., 2012;Colella et al., 2012;Peng & Rubin, 2018). It remains to be seen whether patch-driven models can allow the higher propagation rates seen here and observed by Ghosh et al. (2010), particularly in ruptures that are just a few km wide.
If they cannot, it may be worth considering whether the slow slip rheology varies with time, perhaps because the pore pressure changes (Peng & Rubin, 2017;Rubin, 2011), or whether fault properties vary in space to allow locally reduced resistance to high slip rates.

Potential Consistency With a Slow Earthquake Continuum
It would be particularly interesting to consider spatially variable fault properties if we knew that a single fault zone process produced the entire range of slow earthquakes, from slow slip events to tremor LFEs (Ide et al., 2007). Some rheologies proposed to explain slow slip are unlikely to produce very wide-ranging slip rates, particularly if slip speeds increase as ruptures get smaller. For instance, a rheology where slip rate depends on temperature rather than patch size is unlikely to allow slip rates that increase by a factor of 10,000 as patches get smaller (Hawthorne & Rubin, 2013c;Matsuzawa et al., 2010;Shibazaki & Iio, 2003). Size-limited models, where slip rates tend to be 10-100 times the driving slip rate (Y. J. Liu & Rice, 2005Rubin, 2008;Skarbek et al., 2012;Wei et al., 2018), may also be unlikely to produce very high slip rates.
It is interesting to recall, however, that slip and propagation rates are not just a function of a single, local fault property. The rates may evolve as a function of heterogeneous fault properties and as a function of the current stress field. Indeed, tremor and slow slip propagation sometimes appears diffusive; the propagation slows as a rupture grows (e.g., Amoruso & Crescentini, 2009;Ando et al., 2012;Nakata et al., 2011;Obara, 2012). That diffusive nature could indicate that slow earthquakes start in a region of concentrated stress or on patches that allow high slip rates. They may then spread outward, following Brownian behavior. Standard diffusive models imply that propagation rates scale as V r ∼ T −1/2 , close to the V r ∼ T −2/3 scaling obtained if we assume magnitude-independent stress drops along with a M 0 ∼ T scaling in slow earthquakes, and previous work has reproduced the M 0 ∼ T scaling with diffusion-motivated models (Ide, 2008(Ide, , 2010aIde & Maury, 2018;Ide et al., 2007).
The propagation we observe may provide one more reason to pursue diffusion or patch-driven models, which allow higher slip rates in smaller slow earthquakes. We have added observations of propagation that fill in trends defined by previously observed slow earthquakes ( Figure 1) and thus provide another indication that slow earthquakes constitute a continuum with size-dependent slip rates. The match with previously defined trends is not perfect, but the mismatch could result from observational bias in tremor detection. Our results are restricted by the methodology to a size range between 1 and 8 km, and some others' results are also restricted. For instance, the RTRs and slow slip fronts identified by Houston et al. (2011) and Bletery et al. (2017) are constrained to be longer than 10 km because they used tremor with location spacing or accuracy of 5-10 km.
Our propagation rates are difficult to directly compare with the linear moment-duration scaling found by Ide et al. (2007). However, we can roughly compare the two by plotting black lines in each panel of Figure 1, which assume that (a) slow earthquake moments scale linearly with duration, with a moment rate of 3 × 10 12 N m s −1 , and (b) slow earthquakes have magnitude-independent stress drops Δτ around 30 kPa, as is consistent with a few observations but remains poorly constrained (Bletery et al., 2017;Chestler & Creager, 2017;Hawthorne et al., 2016;Rubin & Armbruster, 2013;Schmidt & Gao, 2010;A. M. Thomas et al., 2018). We assume elliptical ruptures with uniform stress drop and a 3:1 aspect ratio, and we estimate the propagation velocity by dividing the length of the ellipse by the rupture duration. The comparison suggests that our propagation rates also fall roughly along the trend defined by observed slow earthquakes' moments and durations.

Potential Inconsistency With a Slow Earthquake Continuum
However, it is too early to firmly infer that all slow earthquakes are governed by the same fault zone processes. Our results fill one observational gap, but other gaps in the slow earthquake spectrum remain, and some researchers have observed scalings that differ from the overall trend in Figure 1c Farge et al., 2020;Gomberg et al., 2016;Michel et al., 2019;Supino et al., 2020).
Further, one could interpret our observed propagation velocities as evidence against a simple continuum of slow earthquakes. The events we analyze are slightly slower than one would expect after extrapolating the propagation velocities of tremor fronts identified by Bletery et al. (2017) and Houston et al. (2011) (Figure 1b). One could argue that there are slow earthquakes with a wide range of propagation rates and durations, located all over the plot in Figure 1b. Our and others' identified events could simply have sizes that reflect our observational capabilities (Gomberg et al., 2016).
Such observational bias does not seem to explain all the trends in the observed events' sizes. For instance, by extracting durations from the C p time series and propagation from the tremor locations, we should be able to identify at least part of the propagation in longer, faster events, even if they are 20 km across. And if 30-s-long M W 5 earthquakes were common, it would be surprising that they have not yet been spotted. But it may also be surprising, at least to our physical intuition, that a single fault zone process could create slow earthquakes with wide-ranging slip rates, so we must be careful to remember that many events could go unobserved.

Conclusions
We have identified thousands of short bursts of tremor beneath Vancouver Island by employing a phase coherence method developed by Hawthorne and Ampuero (2017) and a set of template LFE waveforms created by Bostock et al. (2012). For 17 bursts, we perform a grid search on the fault plane to track the evolution of tremor and likely slip. We find that these minutes-long events have pulse-like ruptures. They move 1-6 km at speeds of 3 m/s-25 m/s. The events' properties fall roughly, though not quite on, the duration-propagation velocity trend defined by previously observed events. These trends provide further, albeit still inconclusive, evidence that slow earthquakes with a wide range of slip rates are created by the same fault zone processes. In any case, they indicate that any complete physical model of slow slip in Cascadia should reproduce not just events that last weeks, with propagation rates of 0.1 m/s, and subevents that last 3 hr, with propagation rates of 5 m/s, but also subevents that last 2 min, with propagation rates of 20 m/s.

Data Availability Statement
The tremor catalogues created in this study are provided in Supporting Information S1. Zipped text files contain lists of spikes in coherence and the coherence time series. We used seismic data from several networks: the Canadian National Seismograph Network (https://doi.org/10.7914/SN/CN), the POLARIS-BC Network (Nicholson et al., 2005), the Earthscope Plate Boundary Observatory Borehole Seismic Network operated by UNAVCO, the Pacific Northwest Seismic Network (https://doi.org/10.7914/SN/UW), the Earthscope USArray Transportable Array (https://doi.org/10.7914/SN/TA). The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata, and/or derived products used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience and EarthScope (SAGE) Proposal of the National Science Foundation under Cooperative Agreement EAR-1261681. To identify peaks in coherence, we used the find_peaks() function of the Scipy python package (Jones et al., 2001). We collected the data in Figure 1 and Table S3 from a variety of papers, as cited, and from the Slow Earthquake Database at http://www-solid.eps.s.u-tokyo.ac.jp/sloweq/ (Kano et al., 2018).