ARTICLES PUBLISHED ONLINE: 27 APRIL 2015 DOI: 10.1038/NGEO2419 Low friction and fault weakening revealed by rising sensitivity of tremor to tidal stress Heidi Houston* At subduction zones, the level of friction on the deep part of the plate boundary fault controls stress accumulation and release, which govern the transfer of stress to the shallower, locked part of the fault that can slip in a megathrust earthquake. In some subduction zones, the deep fault slips slowly, at speeds far less than shallower, regular earthquakes, and is often accompanied by weak seismic waves called tremor. Tremor and slow slip can be triggered by small stress changes induced by ocean or solid Earth tides. Here I use seismic data combined with calculations of tidal stress to determine the influence of tides on 31,000 tremors generated by six large slow-slip events in Cascadia between 2007 and 2012. I find that the sensitivity of tremor to tidal stresses rises during each slip event, as slip at each spot on the fault accumulates. Specifically, tremor rate is an exponential function of tidal stress, and this exponential sensitivity grows for several days, implying that the fault weakens during slip. I use the relationship between tidal stress and tremor to calculate a coefficient of intrinsic friction for the fault and find values of 0 to 0.1, which indicate that the deep fault is inherently weak. E pisodic tremor and slip (ETS) is a recently discovered phenomenon in which weak seismic signals called tremor accompany slowly migrating slip on a plate boundary interface in slow earthquakes with moment magnitudes up to ∼M7.0 and durations of weeks to months1–3 . Slip in ETS is detectable by GPS and tiltmeters, and features rupture propagation velocities and slip velocities that are, respectively, 4.5 and 6 orders of magnitude smaller than those in regular earthquakes1,4,5 . In comparison with signals from regular earthquakes, seismic radiation from tremor is weak, emergent rather than impulsive, and richer in low frequencies1,6 . ETS tends to rupture plate boundaries below the locked zone that breaks in regular earthquakes7–9 . Furthermore, in several subduction zones10–13 , long-term slow slips with little or no tremor have been observed in and above seismogenic depth ranges, although not in Cascadia. In Cascadia and Japan, where the phenomenon was first identified, ETS occurs unusually regularly7,14 . Slow-slip phenomena range from geodetically detected slip events lasting years to seismically detected individual low-frequency earthquakes (LFEs) lasting tenths of a second3,6 . Their durations seem to be linearly proportional to their seismic moment3 , in contrast to regular earthquakes, whose durations are proportional to the cube root of moment15 . The occurrence and amplitude of tremor and slow slip are modulated by tidal stresses16–21 . Here, I detect, analyse and model changes in tidal sensitivity during ETS to constrain fault weakening and frictional properties. Tremor locations Along the Cascadia subduction zone in Washington and southern Vancouver Island, large ETS events rupture the deep plate interface every 12 to 15 months, propagating with speeds averaging 8 km d−1 (ref. 4). Tremor there has been located by an envelope cross-correlation algorithm22 (see Methods). Tremor locations for all of Cascadia from May 2009 onwards are available from the PNSN Tremor Monitor website. Here I use an augmented catalogue containing ∼31,000 tremor locations from the six major ETS events that ruptured northern Washington and southern Vancouver Island between 2007 and 2012 with magnitudes of 6.5 to 6.8, initiating in January 2007, May 2008, May 2009, August 2010, July 2011 and August 2012. Of the six ETS events, five followed a typical pattern4 , which features downdip initiation near Puget Sound followed by asymmetric propagation, primarily northwestwards (Supplementary Fig. 1). In contrast, in August 2012, ETS initiated updip under mid-Vancouver Island and propagated southeast through the study region. This study examines a change in tidal sensitivity during the several days of slip at each location. I, therefore, accounted for the along-strike ETS propagation with a ‘tremor front’ computed from the specific propagation pattern of each ETS (Methods and Supplementary Fig. 1). I assign each tremor to a group based on when it occurred relative to the arrival of the tremor front. Tremor—and, I infer, slip—typically continues at a given location during an ETS for several days at a decreasing rate. In detail, the behaviour of repeating LFE families suggests that during an ETS, slip, particularly at updip locations, occurs intermittently for several days, and correlates with tidal stressing after an initial uncorrelated burst23,24 . Tremor behind the tremor front occasionally propagates in the opposite direction to the front in rapid tremor reversals (RTRs), which imply weakening of the fault after the main slip front has passed through a region4 . RTRs occur on updip parts of the fault, after the main slip front, and at times of encouraging tidal stress21 . This study determines an average, macroscopic character of the fault for several time periods relative to the inferred arrival of slip, averaged over the entire study region and over six ETS events. Although any specific location or individual ETS may not show clearly the influence of tides, combining data from numerous locations and ETS events overcomes the scatter, reveals a strong influence of tidal stress on tremor occurrence, and even provides a ‘direct image’ of the Coulomb failure condition for frictional sliding. Department of Earth and Space Sciences, University of Washington, Johnson Hall 070, 4000 15th Avenue NE, Seattle, Washington 98195-1310, USA. *e-mail: heidi.houston@gmail.com NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience 1 NATURE GEOSCIENCE DOI: 10.1038/NGEO2419 49.0 Latitude (° N) 49.0 0.70 0.65 48.5 0.60 0.55 48.0 0.50 0.45 47.5 c Initial tremor 0.65 48.5 0.60 0.55 48.0 0.50 0.45 47.5 47.0 0.35 123.5 0.65 48.5 0.60 0.55 48.0 0.50 0.45 47.5 47.0 122.5 Longitude (° W) 2.0 48.5 1.5 48.0 1.0 47.5 0.5 0.40 47.0 0.35 123.5 Later to initial consistency 49.0 0.70 0.40 0.40 d Later tremor 49.0 0.70 Latitude (° N) b All tremor Latitude (° N) a Latitude (° N) ARTICLES 122.5 0.35 123.5 Longitude (° W) 47.0 122.5 Longitude (° W) 123.5 122.5 0.0 Longitude (° W) Figure 1 Consistency of tremor with sense of tidal stress. a, All tremor within the study region. Colour gives fraction of tremor during encouraging tidal stress (Coulomb stress with a frictional coefficient 0.1). Colourbars are centred on 0.53, the fraction of time stress is encouraging. b, Tremor that occurred before 1.5 days after the passage of tremor front. c, Tremor that occurred later than 1.5 days after the passage of tremor front. This group of tremor shows greater consistency with tides. d, Ratio of c to b. Updip regions tend to show greater consistency with tidal stresses and greater increase in consistency. Tidal loading Tidal stresses are generated by gravitational attraction between the Earth and the Sun and Moon. The full stress tensor at depth due to elastic deformations from body tides and ocean loads is estimated20 (see Methods), providing a relatively well-characterized signal with which to assess the effect on earthquakes, tremor or slow slip. Supplementary Fig. 2 shows the grid of locations where tidal stresses are estimated every 10 min in time. Coulomb stress Recent analysis of the effects of tidal stress on tremor on the deep San Andreas Fault concluded that ductile mineral flow mechanisms are inconsistent with the tide–tremor relationship, implying that brittle-frictional behaviour prevails25 . The point of view taken here is that tremor is generated by brittle fracture of small asperities surrounded by frictionally creeping regions. This analysis constrains macroscopic frictional properties of the surrounding creeping regions rather than those of the tremor asperities. Previous studies of the influence of tides on tremor have investigated the effects of shear stress and normal stress separately. For simplicity, here they are combined into Coulomb stress, a physically based combination of stresses that provides a simple way to describe friction: 1CFF = 1τ + µ1σ where 1CFF is the change in Coulomb stress (positive promotes shear failure), 1τ is the change in shear stress on the fault plane in the slip direction, 1σ is the change in normal stress (positive is tensile) and µ is the coefficient of intrinsic friction. A range of friction coefficients is considered, including non-laboratory friction values (that is, low values of µ); weak influence of normal stress implies a low friction coefficient. Pore pressure acts against normal stress and its effect on Coulomb stress should be included. However, it is not known which poroelastic model is most appropriate26,27 . The simplest assumption appropriate for the deep fault interface is that isotropic, undrained conditions prevail. In this case much stress is supported by the pore pressure versus by the rock matrix. Experimental determinations of B in a range of rock types find wide variations, from 0.2 to 0.9 (Table 5 of ref. 28), with, however, few constraints for a deep subduction-zone environment. In such an environment, where ambient pore pressure is high, B = 0.5 is a reasonable choice for the isotropic case. Undrained conditions are consistent with the low permeability thought to prevail on the plate interface29–31 and with the short timescale of tidal loading. With this assumption, my analysis can constrain the intrinsic friction coefficient, because normal and mean stress changes often differ significantly (that is, by amounts greater than the shear stress). An alternative assumption made explicitly or implicitly in some studies of tremor is that highly anisotropic, undrained conditions prevail, in which increases in normal stress are nearly exactly opposed by increases in pore pressure. This assumption has been commonly made in studies of static stress triggering of aftershocks, and was explored by refs 26,27. Then 1CFF = 1τ + µ(1σ + 1p) with 1p = −B1σ (2) 1CFF = 1τ + µ(1 − B)1σ (3) so Such a situation might be approximated if a highly anisotropic fault hosted fluid-filled fractures aligned along the fault zone20 . It is important to note that, in this scenario, the value of B is necessarily close to 1 (refs 20,32). What is constrained in this scenario is the product µ(1 − B), which is effective friction. Uncertainties in B lead to large uncertainties in intrinsic friction µ, hampering the inference of intrinsic friction from effective friction. Furthermore, formulations such as ref. 19 that do not include a poroelastic component are equivalent to equation (3), but with only the product µ(1 − B) constrained. Thus, in ref. 19, the preferred friction value of 0.02 (their Fig. 4) is an effective friction, for which B ranging from 0.9 to 0.96 would imply an intrinsic friction of 0.3 to 0.5. Evolution of tremor sensitivity to tidal stress 1CFF = 1τ + µ(1σ + 1p) with 1p = −B1σm (1) where 1p is the change in pore pressure, 1σm is the change in mean stress and B is Skempton’s coefficient. B is a measure of how 2 The frequency of occurrence and intensity of tremor during ETS are influenced by stresses from solid Earth tides and ocean loading16 . Analysing together all the tremors in the study area demonstrates the known influence of tidal stress on tremor in Cascadia (Fig. 1a). NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience NATURE GEOSCIENCE DOI: 10.1038/NGEO2419 a ARTICLES b 7,000 Number of tremor 5,000 4,000 3,000 2,000 1,000 0 5 Days after tremorfront 10 7,000 15 3.5 3.0 2.5 2.0 1.5 1.0 5.0 −2 µ = 0.1 −1 0 1 2 Coulomb stress change from tides (kPa) 3 B = 0.5 4.5 5,000 Number of tremor 4.0 −3 d 12,188 9,421 9,554 6,000 4,000 3,000 2,000 1,000 0 −5 B = 0.5 0.5 Observed to expected tremor occurrence c µ = 0.1 4.5 Observed to expected tremor occurrence 6,000 0 −5 5.0 12,188 18,975 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0 5 Days after tremorfront 10 −3 15 −2 −1 0 1 2 Coulomb stress change from tides (kPa) 3 Figure 2 Evolution of tremor response to tides. a, Histogram of tremor timing relative to tremor front (note asymmetry). Colours define the grouping and total numbers of tremor. b, Squares plot ratios of observed-to-expected occurrences of tremor versus tidal stress perturbation. Stress bins contain at least 20 tremors and do not overlap. Error bars show standard deviations from bootstrap resampling with replacement. Colourcoded lines show probability distributions of stresses at times of tremor. Exponential function with an exponent 0.351CFF shown for reference. c, As in a but with three groups of tremor. d, As in b but using the tremor grouping in c. Exponential function with an exponent 0.51CFF shown for reference. Sensitivity to tidal stress increases for at least three days as slip accumulates. Of 31,163 tremors, 60% occur during the 53% of time with positive Coulomb stress, for µ = 0.1 and B = 0.5. Assuming each tremor occurrence is independent, the probability that event occurrence is random relative to tidal stressing is negligible, with similar nonrandomness inferred for a wide range of µ and B values. Several studies have noted the influence of tidal stressing on tremor intensity and frequency of occurrence. Tremor activity during three Cascadia ETS fluctuated at periods of 12.5 and 24–25 h, corresponding to the lunar and lunisolar tides16 . For Nankai, Japan and Vancouver Island, Canada, refs 17 and 18, respectively, constructed basic models of tidal stress changes on the plate interface and compared them against the modulation of tremor rates during ETS. The comparisons indicate a tidal influence. On Vancouver Island, tremor was advanced by ∼3 h relative to peak Coulomb stress. Tremor rate at an isolated spot downdip on the Nankai subduction zone correlates exponentially and unusually strongly with low water levels33 . Tremor on the deep San Andreas Fault is strongly influenced by tidal shear stresses19 . Strainmeter signals during Cascadia ETS show that, in addition to modulating tremor, tidal stress also modulates slow slip20 . Tidal stressing affects tremor propagation patterns within ETS as well, inducing rapid tremor reversals during two closely observed Cascadia ETS (ref. 21). Here, for the first time, an evolution of tremor sensitivity to stress perturbations during ETS is seen. Grouping many tremors from different locations and ETS events according to their time of occurrence relative to the tremor front, reveals an increase in sensitivity to tidal stress as slip accumulates over several days at a spot (Figs 1b,c and 2). This indicates a change in the fault condition probably due to weakening. Several groupings of tremor were investigated; two are discussed here. First, tremors are separated into two groups: those before 1.5 days after the tremor front and those after that time (Fig. 2a). This time period of 1.5 days seems optimal, probably because location errors introduce uncertainty in the grouping relative to the tremor NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience 3 NATURE GEOSCIENCE DOI: 10.1038/NGEO2419 ARTICLES a Shear stress change (kPa) Normalized tremor density plotted on stress state 3.0 2 2.5 1 2.0 0 1.5 −1 1.0 −2 −3 −4 −3 −2 −1 0 1 2 3 4 5 (Normal stress + pore pressure) change (kPa) 6 Table 1 Measures of tidal influence on tremor. Tremor group Time after tremor front (d) Sensitivity (kPa−1 ) Consistency Initial −4 to 1.5 0.104 0.548 Middle 1.5 to 3.0 0.152 0.592 0.5 Late 3.0 to 10 0.537 0.679 0.0 Sensitivity—multiplicative constant in the exponential function of stress that best fits ratios of observed-to-expected tremor versus tidal stress (Coulomb stress with µ = 0.1 and B = 0.5). Consistency—fraction of tremor that occurs when tidal stress is encouraging (positive). Should be compared with the fraction of time when stress is positive (0.535). b consistency with the sense of stress (Table 1). Here, the sensitivity of the late group is even greater than that of the late group in Fig. 2a,b, and furthermore, the middle group has a sensitivity intermediate between those of the two groups in Fig. 2a,b. Thus, the sensitivity to stressing evolves over at least ∼3 days (Fig. 2c,d), suggesting that the fault continues to weaken as slip accumulates for several days after the tremor front has passed through a region, even though slip may slow after a fraction of a day23 . Ratio of observed to expected tremor 8 6 4 2 0 −6 −4 −2 0 2 4 6 1 2 3 8 −4 −3 −2 −1 0 Shear stress change (kPa) 4 Figure 3 Influence on tremor of shear versus normal stress changes. a, Ratio of observed-to-expected tremors (late group) plotted in effective-normal-stress versus shear-stress space. Effective normal stress is normal stress plus pore pressure from Model 1 with B = 0.5 (see equation (1)). Expected number of tremors is based on the time the fault occupies that stress state. Stress bins do not overlap. Positive normal stress is tensile. Dashed red line has slope −0.1, indicating a very shallow slope of the Coulomb sliding line—that is, low intrinsic friction. b, 3D perspective view of a, showing its exponential character. front of the slowly moving ETS pulse. The two groups of tremor exhibit very different responses to the tidal stresses, as demonstrated by normalized histograms (Fig. 2b). They show the fraction of tremor that occurred in a given Coulomb stress range divided by the fraction of time the stress at the nearest gridpoint occupied that stress range, for several days around the tremor (Methods). Here Coulomb stress was calculated with a low µ, as justified below, but a similar distinction between the early and later groups is seen at levels of µ up to 0.6. Figure 1 shows the spatial variations in the ‘consistency’, the fraction of tremor consistent with positive tidal stresses. For the initial tremor group, consistency is relatively high in the region where of five of the six ETS events initiate (Fig. 1b and Supplementary Fig. 1), which suggests that before ETS the stress there is closer to failure than in other parts of the region. As noted above, tidal influence is much greater on the second, later group (Fig. 1c). In general, it is greater in updip regions, and near major waterways and Vancouver Island (see Supplementary Fig. 3). Figure 1d shows the change in consistency between the two groups (that is, the fractions in Fig. 1c divided by those in Fig. 1b). Locations more consistent with, and therefore more responsive to, tidal stressing could mark zones with a different rheology or particularly high pore pressure. Locations with larger changes after the initial slip pulse may mark sites of greater weakening. Second, tremors are divided into three groups to further analyse the evolution of sensitivity over time (Fig. 2c,d). For late tremor, the marked increase in the ratio strongly suggests an exponential function of stress. The exponent factor is a measure of tidal response, which I call the ‘sensitivity’, in addition to the above-defined 4 Direct image of low-friction Coulomb sliding line A closer look at the relationship between tremor and stress state is afforded by plotting the normalized frequency of occurrence of tremor (in this case, for the last and most sensitive group) on a base of effective normal versus shear stress (Fig. 3). Here again, the frequency of occurrence is normalized for each stress state by the amount of time the fault occupies that stress state over the entire several days when tremor is in the vicinity. Remarkably, a shallow-sloping feature related to a Coulomb frictional sliding failure line is directly imaged in Fig. 3. The slope implies an intrinsic friction of 0.1, much smaller than laboratory-derived estimates of approximately 0.6 (that is, Byerlee friction). The exponential rise of the normalized frequency of tremor with increasing Coulomb stress can be seen in the perspective view in Fig. 3b. Thus, Fig. 2d can be visualized as a stack along the Coulomb frictional sliding line. Implications for intrinsic friction Whether the tremor–tide correlation has the potential to constrain intrinsic friction depends on which poroelastic model is more applicable. The formulation for Coulomb stress using equation (1) (termed Model 1) retains intrinsic friction, whereas in equation (2) (termed Model 2) intrinsic friction is multiplied by (1 − B), where B is Skempton’s coefficient. Furthermore, if Model 2 applies, then the value of B must be close to 1, as discussed next, and its precise value is not well known. Thus, in the formulation of Coulomb stress with Model 2, intrinsic friction trades off with B and typically cannot be constrained. Poroelastic Model 2 implies a Skempton’s coefficient B of close to 1. Model 2 could describe an extremely anisotropic fabric, in which fluid-filled cracks enlarge in only one direction—normal to the fault plane20,32 . Such behaviour is approached far from the edges of large, flat cracks with a small ratio of opening-width to radius (a ratio of 1 to 100 suggested in ref. 20). Such cracks are proposed to be aligned with the fault plane and might develop in association with the large displacements and possible geometrical and petrological complexity in the fault zone. Model 2 describes the poroelastic response for this extreme endmember. The Skempton’s coefficient for Model 2 must be close to 1 (as in ref. 20 and Fig. 2 of ref. 32) because a fabric that produces strain only normal to the fault plane, such as large flat cracks covering over 90% of the fault plane20 , also supplies very little resistance of the solid matrix to normal opening under fluid pressure. NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience NATURE GEOSCIENCE DOI: 10.1038/NGEO2419 ARTICLES a b B = 0.4 B = 0.5 B = 0.6 B = 0.1−0.9 40 10 38 36 0 Percentage excess tremor Percentage excess tremor 5 More tremors than expected Fewer tremors than expected −5 −10 34 32 30 28 26 24 22 −15 0.4 0.5 0.6 0.7 0.8 0.9 Skempton’s coefficient, B 20 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Intrinsic friction coefficient Figure 4 Pore-pressure models and friction coefficient. a, Excess observed-to-expected tremor for the subset of times most sensitive to difference between poroelastic Models 1 (purple) and 2 (green). For each Skempton’s coefficient B, µ-values from 0.1 to 0.6 were tested. Larger triangles indicate probability <0.01 of random occurrence from χ 2 analysis. Error bars are standard deviations. Most results for Model 2 give fewer tremors than expected in the absence of tidal influence. Purple box encloses best results for Model 1 (preferred), green box encloses best results for Model 2. b, Using Model 1 with B = 0.5 and the full data set of late-group tremor above a minimum stress of 0.35 kPa, the tremor–tide correlation is maximized for intrinsic friction from 0 to 0.1. However, tremor correlates better with tidal stresses computed with Model 1, in which pore pressure varies with mean, rather than normal, stress (Fig. 4a). The correlation of the entire tremor data set with tidal stresses is only slightly stronger for Model 1 than Model 2. But, if one considers a subset of the data sensitive to the differences between the models, then a significant distinction emerges. Such a subset consists of the times when Models 1 and 2 yield opposite signs for the Coulomb stress change. Shear stress tends to be small for this subset of times (so that the two models can change the sign of the Coulomb stress), and there are relatively few tremors, particularly for B < 0.4. The percentage excess tremor of this sensitive subset was computed for the two models and a range of friction and Skempton’s coefficients (Fig. 4a). The percentage excess for Model 1 is positive, and mostly higher than that for Model 2, which yields a negative tremor–tide correlation for most cases. Preferred Skempton’s coefficients for the two models are outlined by boxes in Fig. 4a. Thus, Model 1 produces a stronger tremor–tide correlation, and therefore probably better represents the actual poroelastic behaviour than Model 2, an extremely anisotropic endmember. Therefore, the tremor–tide correlation can yield constraints on intrinsic friction. Such constraints are best derived by analysis of the full data set, as in the previous sections. Figure 4b implies that the coefficient of intrinsic rock friction deep on the Cascadia plate interface is less than 0.15—far smaller than laboratoryderived estimates (that is, Byerlee friction). This result depends on the adoption of Model 1, not on the specific value assumed for B (Fig. 4b). Several instances of very low friction on a subductionzone interface have been recently reported34 . Frictionally weak minerals such as talc35–37 , saponite30 , or graphite38 may play a role in controlling fault weakening and creep in mature largedisplacement fault zones. Fault gouge from the 3-km-deep drill hole on the SAF contains thin authigenic hydrous clay coatings that may facilitate creep, which occurs when the coatings are sufficiently interconnected39 . Modelling the evolution of stress and strength during ETS A stress threshold failure model gives insight into the increasing influence of tidal stress on tremor during slip at a spot, as well as the exponential behaviour at later times (seen in Fig. 2b,d). Rather than distinguishing between creeping regions and brittle asperities, the model posits a macroscopic fault region that encompasses both, creeping frictionally when stress exceeds strength, generating tremors immediately. The strengths are randomly drawn from a Weibull distribution every ten minutes for five days. Stress on the fault is modelled as the sum of slow tectonic loading, estimated transient stress from the propagating ETS slip pulse, and the actual oscillating tidal Coulomb stress at that location (Fig. 5a). If the stress exceeds the strength, a tremor occurred. Parameters are chosen based on various observational constraints (see Methods). Ambient or background stress levels following ETS are not known and not needed for this model. Figures 2b,d and 3 imply that they are at least greater than the tidal fluctuations of a few kPa, as the larger negative shear or Coulomb tidal stresses are not sufficient to trigger tremor. This threshold failure model matches the relationship between tremor and stress amplitude well (compare Figs 2d and 5b), and also explains the first-order features of observed time shifts between tremor and peak stress (see Supplementary Information). The observed exponential behaviour is generated by the model because stress levels late in the five-day period are oscillating over the lower half of a Weibull distribution, which falls off exponentially. The mean strength must fall during ETS by roughly the amount that the stress drops, because tidal oscillations are small relative to that drop (Fig. 5a). Such a strength decrease might result from breakage of mineral precipitates40 built up over the ∼14 month recurrence interval. The timing and amount of weakening inferred here provide new constraints for modelling of slow slip and tremor. The low level of macroscopic intrinsic friction inferred is the first in situ measurement of intrinsic friction on the deep plate interface. NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience 5 NATURE GEOSCIENCE DOI: 10.1038/NGEO2419 ARTICLES a 120 Slip at a spot during ETS Stress or strength (kPa) 100 80 Stress 60 Strength distribution 40 20 0 9 10 11 12 13 14 15 16 17 18 Time (days) b Ratio of simulated to expected tremor 5 4 3 2 1 0 −3 −2 −1 0 1 2 3 Coulomb stress change from tides (kPa) Figure 5 Model of stress and strength evolution during ETS. a, Solid line shows stress history at a spot for one realization of threshold failure model. Stress is the sum of tidal stresses, slow tectonic loading, and transient stress from slip pulse. Colour indicates time periods relative to tremor front, as in Fig. 2c. Every 10 min, strength is drawn from a Weibull distribution. Circle highlights where the exponential part of the distribution interacts with oscillating tidal stress, enhancing tremor exponentially. Absolute stress levels are arbitrary. b, Tremor occurrence versus tidal stress from 500 model realizations with randomly chosen locations and five-day intervals. Note the similarity to Fig. 2d. Methods Methods and any associated references are available in the online version of the paper. Received 29 December 2014; accepted 18 March 2015; published online 27 April 2015 References 1. Obara, K. Nonvolcanic deep tremor associated with subduction in southwest Japan. Science 296, 1679–1681 (2002). 2. Rogers, G. & Dragert, H. Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip. Science 300, 1942–1943 (2003). 3. Ide, S., Beroza, G., Shelly, D. & Uchide, T. A scaling law for slow earthquakes. Nature 447, 76–79 (2007). 6 4. Houston, H., Delbridge, B., Wech, A. & Creager, K. Rapid tremor reversals in Cascadia generated by a weakened plate interface. Nature Geosci. 4, 404–409 (2011). 5. Rubin, A. Designer friction laws for bimodal slow slip propagation speeds. Geochem. Geophys. Geosyst. 12, Q04007 (2011). 6. Shelly, D., Beroza, G. & Ide, S. Non-volcanic tremor and low-frequency earthquake swarms. Nature 446, 305–307 (2007). 7. Hirose, H. & Obara, K. Short-term slow slip and correlated tremor episodes in the Tokai region, central Japan. Geophys. Res. Lett. 33, L17311 (2006). 8. Dragert, H., Wang, K. L. & James, T. S. A silent slip event on the deeper Cascadia subduction interface. Science 292, 1525–1528 (2001). 9. Nishimura, T., Matsuzawa, T. & Obara, K. Detection of short-term slow slip events along the Nankai Trough, southwest Japan, using GNSS data. J. Geophys. Res. 118, 3112–3125 (2013). 10. Hirose, H. & Obara, K. Repeating short- and long-term slow slip events with deep tremor activity around the Bungo Channel region, southwest Japan. Earth Planets Space 57, 961–972 (2005). 11. Miyazaki, S., Segall, P., McGuire, J. J., Kato, T. & Hatanaka, Y. Spatial and temporal evolution of stress and slip rate during the 2000 Tokai slow earthquake. J. Geophys. Res. 111, B03409 (2006). 12. McCaffrey, R., Wallace, L. M. & Beavan, J. Slow slip and frictional transition at low temperature at the Hikurangi subduction zone. Nature Geosci. 1, 316–320 (2008). 13. Ozawa, S., Suito, H. & Tobita, M. Occurrence of quasi-periodic slow-slip off the east coast of the Boso Peninsula, Central Japan. Earth Planets Space 59, 1241–1245 (2007). 14. Brudzinski, M. R. & Allen, R. M. Segmentation in episodic tremor and slip all along Cascadia. Geology 35, 907–910 (2007). 15. Houston, H. Influence of depth, focal mechanism, and tectonic setting on the shape and duration of earthquake source time functions. J. Geophys. Res. 106, 11137–11150 (2001). 16. Rubinstein, J., La Rocca, M., Vidale, J., Creager, K. & Wech, A. Tidal modulation of nonvolcanic tremor. Science 319, 186–189 (2008). 17. Nakata, R., Suda, N. & Tsuruoka, H. Non-volcanic tremor resulting from the combined effect of Earth tides and slow slip events. Nature Geosci. 1, 676–678 (2008). 18. Lambert, A., Kao, H., Rogers, G. & Courtier, N. Correlation of tremor activity with tidal stress in the northern Cascadia subduction zone. J. Geophys. Res. 114, B00A08 (2009). 19. Thomas, A., Nadeau, R. & Burgmann, R. Tremor-tide correlations and near-lithostatic pore pressure on the deep San Andreas Fault. Nature 462, 1048–1051 (2009). 20. Hawthorne, J. & Rubin, A. Tidal modulation of slow slip in Cascadia. J. Geophys. Res. 115, B007502 (2010). 21. Thomas, T. et al. Evidence for tidal triggering of high-amplitude rapid tremor reversals and tremor streaks in northern Cascadia. Geophys. Res. Lett. 40, 4254–4259 (2013). 22. Wech, A. G. & Creager, K. C. Automated detection and location of Cascadia tremor. Geophys. Res. Lett. 35, L20302 (2008). 23. Sweet, J. R. Unlocking the secrets of slow slip in Cascadia using low-frequency earthquakes. PhD thesis, Univ. Washington (2014). 24. Royer, A. A., Thomas, A. M. & Bostock, M. G. Tidal modulation and triggering of low-frequency earthquakes in northern Cascadia. J. Geophys. Res. 120, 384–405 (2015). 25. Beeler, N. M., Thomas, A., Burgmann, R. & Shelly, D. Inferring fault rheology from low-frequency earthquakes on the San Andreas. J. Geophys. Res. 118, 5976–5990 (2013). 26. Beeler, N. M., Simpson, R. W., Hickman, S. H. & Lockner, D. A. Pore fluid pressure, apparent friction, and Coulomb failure. J. Geophys. Res. 105, 25533–25542 (2000). 27. Cocco, M. & Rice, J. Pore pressure and poroelasticity effects in Coulomb stress analysis of earthquake interactions. J. Geophys. Res. 107, 2030 (2002). 28. Paterson, M. S. & Wong, T-f. Experimental Rock Deformation—The Brittle Field 2nd edn (Springer, 2005). 29. Ingebritsen, S. E. & Manning, C. E. Permeability of the continental crust: Dynamic variations inferred from seismicity and metamorphism. Geofluids 10, 193–205 (2010). 30. Sone, H., Shimamoto, T. & Moore, D. E. Frictional properties of saponite-rich gouge from a serpentinite-bearing fault zone along the Gokasho-Arashima Tectonic Line, central Japan. J. Struct. Geol. 38, 172–182 (2012). 31. Audet, P., Bostock, M. G., Christensen, N. I. & Peacock, S. M. Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing. Nature 457, 76–78 (2009). 32. Endres, A. L. Geometrical models for poroelastic behaviour. Geophys. J. Int. 128, 522–532 (1997). 33. Ide, S. & Tanaka, Y. Controls on plate motion by oscillating tidal stress: Evidence from deep tremors in western Japan. Geophys. Res. Lett. 41, 3842–3850 (2014). NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience NATURE GEOSCIENCE DOI: 10.1038/NGEO2419 34. Fulton, P. M. et al. Low coseismic friction on the Tohoku-Oki Fault determined from temperature measurements. Science 342, 1214–1217 (2013). 35. Moore, D. E. & Rymer, M. J. Talc-bearing serpentinite and the creeping section of the San Andreas Fault. Nature 448, 795–797 (2007). 36. Moore, D. E. & Rymer, M. J. Correlation of clayey gouge in a surface exposure of serpentinite in the San Andreas Fault with gouge from the San Andreas Fault Observatory at Depth (SAFOD). J. Struct. Geol. 38, 51–60 (2012). 37. Hirauchi, K., den Hartog, S. A. M. & Spiers, C. J. Weakening of the slab-mantle wedge interface induced by metasomatic growth of talc. Geology 41, 75–78 (2013). 38. Oohashi, K., Hirose, T. & Shimamoto, T. Graphite as a lubricating agent in fault zones: An insight from low- to high-velocity friction experiments on a mixed graphite-quartz gouge. J. Geophys. Res. 118, 2067–2084 (2013). 39. Schleicher, A. M., van der Pluijm, B. A. & Warr, L. N. Nanocoatings of clay and creep of the San Andreas Fault at Parkfield, California. Geology 38, 667–670 (2010). ARTICLES 40. Audet, P. & Burgmann, R. Possible control of subduction zone slow-earthquake periodicity by silica enrichment. Nature 510, 389–392 (2014). Acknowledgements I thank J. Hawthorne for assistance with calculations of tidal stresses, S. Ide for suggesting the use of a Weibull distribution, and J. Vidale and K. Creager for helpful suggestions. Reviews from R. Burgmann and the Editor helped improve the paper. This work was supported by NSF grant EAR-1447005. 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 author declares no competing financial interests. NATURE GEOSCIENCE ADVANCE ONLINE PUBLICATION www.nature.com/naturegeoscience 7 ARTICLES Methods Tremor locations. Tremors were located by an envelope cross-correlation algorithm applied to overlapping 5-min long windows of envelopes of seismograms22 . Epicentral location errors are estimated at ∼8 km. Depth is poorly constrained and the tremors are inferred to lie on a plate interface model41 . Most available evidence supports the location of tremors very near to or on the plate interface, but the exact relationship of tremors to plate interface remains disputed42–44 . Calculation of the tremor front. Tremor epicentres are rotated to slab coordinates which are defined relative to a reference line fit to the updip edge of tremor along the Cascadia subduction zone22 . As ETS propagation is mainly along strike, with lengths of several hundred km and widths of ∼50 km, epicentres are grouped into overlapping 4-km-long bins in the along-strike direction. The tremor front is then defined as the time when 5% of the tremor that ultimately occurs in a bin has already occurred, leaving 95% still to come (Supplementary Fig. 1). Before that step, a correction is applied to account for systematic timing differences due to propagation over the length of the bin. This approach allows easy handling of bilateral propagation. Tidal stress calculations. Stresses from elastic deformations due to both water loads and bodily tides in the solid Earth are calculated following ref. 20, which uses the software package SPOTL (ref. 45) and integrates strain at depth from water loads with spatially varying height, including the Pacific Ocean and regional waterways (including Puget Sound). Water loads within five degrees of the calculation location are included. In contrast, the solid Earth body tides have a much longer spatial wavelength, permitting the use of strains computed at Earth’s surface. Both types of strains are converted to stresses using a Poisson’s ratio of 0.25 and an elastic modulus of 30 GPa, and summed. Time-varying stress tensor components are computed on a curving spatial grid on the inferred plate interface at depth41 , with ∼12 km spacing (Supplementary Fig. 2). The resulting full stress tensors are projected to normal and shear (in the plate convergence direction) stresses on the interface, as well as mean stress (one-third the tensor trace). Five tidal periods are included: M2 , S2 , N2 , K1 and O1 . Time series of these stresses are computed every 10 min for the duration of each ETS. Normalization. Histograms are constructed with ratios of numbers of observed-to-expected tremors in each stress bin. The number of tremors that would be expected in the absence of sensitivity to tidal stresses is calculated separately for each grid element to avoid artefacts due to spatially varying tremor densities and total stress ranges. For a given stress bin, the expected number of tremors in a grid element is the fraction of time that Coulomb stress values there lie within that stress bin times the total number of tremors observed there. Using all stresses in that grid element during a several-day-interval around tremor to estimate expected tremors, rather than several years of stresses, improves accuracy, because the distribution of stresses at non-ETS times is not relevant here. Coulomb stress values are computed with shear stress calculated in the convergence direction on the plate interface. For the analysis in Fig. 4b, specifying a positive minimum stress level includes only larger, better-resolved stresses and achieves a more robust result. Threshold failure model and parameters. For a given fault gridpoint, a random 5-day period of actual oscillating tidal stresses at that location is chosen. NATURE GEOSCIENCE DOI: 10.1038/NGEO2419 Every 10 min, the tidal, transient and tectonic stresses are summed and compared against randomly drawn strengths, with a tremor generated if stress exceeds strength. This is repeated 500 times at random fault gridpoints with random 5-day intervals of actual tidal stresses. Model parameters, although non-unique, are constrained by various considerations. Geodetic inversions of major ETS in this region typically find 2–3 cm of slip distributed over ∼50 km downdip width, resulting in stress drops of 15–30 kPa (refs 46,47). The tectonic loading rate is estimated assuming steady stress recharge by plate convergence. It is orders of magnitude below the tidal or transient loading rates. The amplitude of stress from the transient slip pulse is roughly based on inferences from elasticity and kinematics that the peak is several times the subsequent stress drop5 . The duration of the transient is based on the tide–tremor observations presented above (Fig. 2d). Heterogeneity of material breaking strength is commonly characterized by a Weibull distribution with a mean, scale parameter and shape parameter48 . Parameter values are based on those determined for material failure (that is, shape parameter between two and six48 ). The mean of the strength distribution must decrease to allow tremor during, but not before, the transient slip pulse, as well as during the last three days of slip at a location. In Fig. 5a, the mean of the strength distribution decreases by 15 kPa over ∼3 days, the scale and shape parameters are 10 kPa and 3, respectively. The resulting tremors are divided into three groups relative to the start of the transient at each gridpoint (that is, the start of the 5-day period, close to the tremor front) and plotted in Fig. 5b as in Fig. 2d. The similarity of the tidal sensitivities is notable. Code availability. Codes used to calculate tidal stress are available at http://igppweb.ucsd.edu/~agnew/Spotl/spotlmain.html. We have opted not to make the other codes associated with this work available owing to their specificity and complexity. References 41. McCrory, P., Blair, J., Waldhauser, F. & Oppenheimer, D. Juan de Fuca slab geometry and its relation to Wadati-Benioff zone seismicity. J. Geophys. Res. 117, B009407 (2012). 42. Shelly, D., Beroza, G., Ide, S. & Nakamula, S. Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip. Nature 442, 188–191 (2006). 43. La Rocca, M. et al. Cascadia tremor located near plate interface constrained by S minus P wave times. Science 323, 620–623 (2009). 44. Kao, H. et al. A wide depth distribution of seismic tremors along the northern Cascadia margin. Nature 436, 841–844 (2005). 45. Agnew, D. C. NLOADF: A program for computing ocean-tide loading. J. Geophys. Res. 102, 5109–5110 (1997). 46. Schmidt, D. A. & Gao, H. Source parameters and time-dependent slip distributions of slow slip events on the Cascadia subduction zone from 1998 to 2008. J. Geophys. Res. 115, B00A18 (2010). 47. Wech, A. G., Creager, K. C. & Melbourne, T. I. Seismic and geodetic constraints on Cascadia slow slip. J. Geophys. Res. 114, B10316 (2009). 48. Wong, T-f., Wong, R. H. C., Chau, K. T. & Tang, C. A. Microcrack statistics, Weibull distribution and micromechanical modeling of compressive failure in rock. Mech. Mater. 38, 664–681 (2006). NATURE GEOSCIENCE www.nature.com/naturegeoscience