## Abstract

The relative motion of tectonic plates is accommodated at boundary faults through slow and fast ruptures that encompass a wide range of source properties. Near the Parkfield segment of the San Andreas fault, low-frequency earthquakes and slow-slip events take place deeper than most seismicity, at temperature conditions typically associated with stable sliding. However, laboratory experiments indicate that the strength of granitic gouge decreases with increasing temperature above 350°C, providing a possible mechanism for weakening if temperature is to vary dynamically. Here, we argue that recurring low-frequency earthquakes and slow-slip transients at these depths may arise because of shear heating and the temperature dependence of frictional resistance. Recurring thermal instabilities can explain the recurrence pattern of the mid-crustal low-frequency earthquakes and their correlative slip distribution. Shear heating associated with slow slip is sufficient to generate pseudotachylyte veins in host rocks even when fault slip is dominantly aseismic.

## INTRODUCTION

Recent seismo-geodetic observations have illuminated episodic deformation at the root of active strike-slip faults and megathrusts alike around the world, revealing a wide variety of slip events that fill the spectrum between slow and fast slip (*1*). Tremor emissions below the seismogenic zone are often associated with slow slip (*2*), but the underlying physics remains elusive. Several mechanisms have been proposed to explain the slow-slip phenomenon, including large nucleation size (*3*), dilatant strengthening (*4*), and semibrittle creep (*5*), but the important role of fluids is often invoked (*6*). The simultaneous tremor or low-frequency earthquake (LFE) emissions are thought to represent a form of deterministic chaos emerging from the nonlinear dynamics (*7*) or the response of a multiscale rock assembly (*8*).

Swarms of LFE can be found on the San Andreas fault near the Parkfield segment (*9*), occupying a separate depth interval in the crust from regular seismicity (Fig. 1), presenting along-strike variations in recurrence patterns, amplitude, and sensitivity to tidal stress (*10*). The shallowest source of mid-crustal LFE, northwest of Parkfield, that recurs frequently, every 1 to 15 months, clustering seismic events for days to weeks at a time, has been clearly associated with slow-slip transients based on geodetic data (*11*). These observations challenge our fundamental assumptions about the rheology of the continental crust. Although the down-dip segmentation of rupture style is widespread at active faults (*12*), these seismic emissions take place at depths typically associated with stable creep in a continental setting, preventing stress accumulation and the development of frictional instabilities.

Rapid slip below the seismogenic zone has been found in different seismic settings (*13*, *14*), often explained in a top-down model where seismic activity in the upper crust drives afterslip and aftershocks at greater depth. In contrast, we argue that the coupling between shear heating and the temperature dependence of frictional resistance and contact healing may be responsible for the development of slow-slip events in the velocity-strengthening domain below the seismogenic zone. Friction laboratory experiments on granitic rocks and quartz gouge under high pore-fluid pressure (*15*, *16*) indicate velocity-weakening and temperature-strengthening friction between about 50° and 350°C, explaining the location of the seismogenic zone (*17*). The same set of experiments also indicates simultaneous velocity-strengthening and temperature-weakening behavior below the seismogenic zone, at ambient temperatures above 350°C. As temperature may vary markedly with shear heating (*18*, *19*), the temperature-weakening behavior of steady-state friction may offer the conditions for episodic slip.

To support these claims, we first describe a constitutive framework for rate-, state-, and temperature-dependent friction, highlighting the necessary condition for the emergence of thermal instabilities with velocity-strengthening behavior. We then describe the emergence of slow-slip events in a simplified spring-slider system in nonisothermal condition. Last, we explain the recurrence pattern and the correlative slip distribution of the shallow LFE sources using three-dimensional simulations. Once tuned to seismo-geodetic observations, the model allows us to discuss the dynamic range of temperature during the seismic cycle and its implications for fault fabric in the middle crust, below the seismogenic zone.

## RESULTS

### Thermally activated constitutive law for fault slip

We simulate the coupled evolution of slip and temperature along the San Andreas fault using a microphysical model of rate-, state-, and temperature-dependent friction based on the evolution of the real area of contact (*20*), whereby the frictional resistance takes the form

(1)where τ and

$\overline{\mathrm{\sigma}}$ are the shear and effective normal stress on the fault; *V* and θ are slip velocity and age of asperity contact, respectively; μ_{0}, *V*_{0}, and *T*_{0} are reference friction coefficient, velocity, and temperature, respectively; *a* ≪ 1 and *b* ≪ 1 are power exponents; and *Q* is an activation energy for the direct effect of temperature. The age of contact is treated as a state variable with the thermally activated evolution law

(2)where *L* is a characteristic weakening distance over which the system evolves toward steady state and *H* is the activation enthalpy for contact healing. An increase in temperature leads to an immediate reduction in frictional strength called the direct effect, but also an acceleration of healing and time-dependent strengthening, leading to a competition between the two opposing effects on steady-state strength (*15*). In isothermal condition with *T* = *T*_{0}, Eqs. 1 and 2 reduce to the multiplicative form of rate-and-state friction (*20*) with the aging law (*21*).

Temperature is allowed to vary dynamically because of shear heating and heat diffusion. We consider that deformation occurs in an active shear zone with a ~1-m-thick fault core made of several subparallel, partially overlapping primary slip surfaces that localize slip over a ∼1-mm-thick active shear layer. The localization of fault slip during any single slip event makes shear heating efficient, even though the overall shear zone may be quite large when considering all the primary and secondary slip surfaces left behind from previous events. As the primary and secondary slip surfaces may overlap partially, the temperature gradient is lower within the fault core than outside of it, and temperature diffusion within the fault core is inefficient. In contrast, because of their high fracture density, the damaged rocks outside the fault core often have a higher fault-perpendicular permeability, and they may advect heat away from the fault effectively (*22*). In addition, because of the widely different length scales involved, thermal diffusion in the fault-parallel direction can be neglected. To keep the problem tractable, we may approximate this structural setting by considering a membrane diffusion where fault temperature can be treated as another state variable with the evolution law

(3)where *T* is the temperature in the shear zone, *D* is the average thermal diffusivity in the surrounding damage zone of thickness *W*, the product ρ*c* is the specific heat per unit volume in the fault zone, *T*_{b} is the bath temperature outside the fault zone, which is considered in drained, isothermal condition, and *w* is the thickness of the active shear layer (Fig. 2). We neglect the evolution of pore pressure on the fault and its effect on strength. Equations 1 to 3 form a complete constitutive framework for the evolution of friction and fault slip. The dynamics of the system is eventually controlled by the spatial distribution of the thermomechanical properties and the mechanisms of stress distribution (*7*).

### Slow-slip events as thermal instabilities

We first explore the range of physical conditions for the emergence of slow-slip events with duration of weeks and recurrence times of the order of 1 year using a spring-slider assembly controlled by Eqs. 1 to 3. The emergence of frictional instabilities is conditioned on the rate and temperature dependence of friction at steady state (*21*, *15*). Velocity weakening occurs for *a* − *b* < 0 and temperature strengthening for *aQ* − *bH* < 0. To explain the seismic emission below the seismogenic zone consistent with laboratory data on granitic rocks and quartz gouge, we focus on a velocity-strengthening (*a* − *b* > 0) but temperature-weakening (*aQ* − *bH* > 0) fault patch embedded in a rock matrix at temperature *T*_{b} = 400°C, representative of mid-crustal conditions with a 20°C/km geothermal gradient. Slow-slip events emerge for a wide range of thermomechanical parameters, but we present a case where the peak temperature throughout the seismic cycle remains below the liquidus of wet granite, between 900° and 1100°C depending on composition. The physical parameters are summarized in table S1.

Away from initial conditions, when the system is continuously loaded at a constant rate, the mechanical system oscillates in closed orbits formed by any pair among velocity, friction, age of contact, and temperature, in repeating cycles (Fig. 2). The system traverses six stages with distinct velocity, temperature, or stress patterns. A long nucleation initiates at stage 1, when stress accumulates, but diffusion continues to dominate the temperature evolution. At the onset of stage 2, shear stress is sufficient for shear heating to overcome diffusion and for temperature to start increasing. At stage 3, temperature weakening accelerates fault slip, leading to a peak in slip velocity when steady state is reached. At this point, temperature is high enough to make healing dominate the state evolution, which strengthens the contact and reduces slip velocity. A peak in temperature is reached at the end of stage 4. At the onset of stage 5, shear stress is low enough for diffusion to dominate the temperature evolution. The transition between stages 4 and 5 is a temperature inflection point. During stage 6, the temperature returns to background levels, the fault relocks, and the stress starts slowly accruing again. As velocity and temperature covary, temperature is often out of phase, with, for example, delayed cooling relative to slip deceleration. These stages unfold in a strictly periodic manner because of the absence of any rheological contrast. The mechanical response of a spring-slider assembly with temperature-weakening, velocity-strengthening friction contrasts with the case of isothermal, velocity-weakening friction as the accelerated healing and restrengthening associated with high temperature in the former allows repeat cycles of slow-slip events without oversized radiation damping.

This simple model demonstrates how the influence of temperature on contact healing and the strength of frictional sliding may govern the nucleation, propagation, and arrest of deep slow-slip events. More complexity in recurrence times, moment, and duration emerges from more sophisticated, three-dimensional simulations.

### Deep slow-slip events on the San Andreas fault

Building upon insights from a simple model, we now simulate the deep slow-slip events along the San Andreas fault on a planar, finite fault embedded in an elastic half-space using the boundary integral method (*7*, *23*). We use the radiation-damping approximation (*24*), as neglecting the wave-mediated stress transfer is appropriate in the slow-slip regime.

We consider a rectangular domain from 10- to 26-km depth extending from 12 to 58 km northwest of Parkfield with velocity-strengthening, temperature-weakening properties throughout (Fig. 3). The thermomechanical properties are uniform, except for the thickness of the active shear layer and the background temperature. We assign a narrow active shear zone thickness of 0.14 mm in a rectangular patch centered around the shallowest LFE sources to promote unstable behavior and the development of thermal instabilities in this region (Fig. 3). The adopted shear layer thickness is on the low range of field observations (*25*). Outside the unstable region, we use a much larger value to promote stability. For the background temperature, we use a thermal gradient of 21°C/km with a surface temperature of 15°C. The physical parameters are summarized in table S2. A low effective normal stress

= 20 MPa is required to reproduce the duration and recurrence times of the LFE bursts, indicating near-lithostatic pore-fluid pressure. In general, low normal stress hinders efficient shear heating, but this is compensated for in the narrow unstable region by a thin active shear zone. We simulate the dynamics of fault slip on the deep San Andreas fault for a period of 300 years and focus on a representative decade that showcases a wide range of rupture sizes. The long simulation time allows us to examine the long-term space-time clustering of slow-slip events and to mitigate the possible bias from initial conditions.

We obtain a complex sequence of slow-slip events associated with spontaneously emerging thermal instabilities. Large variations of cumulative slip distributions can be found among slow-slip events, as rupture may sometimes propagate far into the surrounding stable region because of the smooth transition from slow slip to afterslip and then relocking near the source region (Fig. 3). Simulation of surface displacements at Global Positioning System (GPS) stations during the seismic cycle indicates a weak contribution from the slow-slip events, accounting to less than 1 mm in a decade, explaining the challenge posed by their geodetic detection (*11*). The coupled evolution of slip, stress, and temperature during the seismic cycle (Fig. 4) shows some remarkable differences with the simpler spring-slider model, as more complexity in recurrence patterns can take place in a three-dimensional model. The sequence of slow-slip events is aperiodic with full and partial ruptures of varying sizes of the unstable patch. The slip per event averages 30 μm distributed over 3 km, leading to stress drops of the order of 1 kPa, small enough for the slow-slip cycle to be perturbed by tidal stresses, which are of the order of 0.1 and 4 kPa for the shear and normal components, respectively (*26*).

The numerical simulation generates a large catalog of moment magnitude (Mw) 2.4 to 4.72 events, with a concentration of events with Mw ∼4.5 (Fig. 5), in accordance with the average magnitude of geodetically determined events (*11*). However, the model suggests a wide range of event sizes because of the elongated shape of the unstable region. The simulated slow-slip events recur every few weeks to about 20 months, most frequently between 1 and 4 months, with a duration between 1 and 5 weeks, with most events propagating for 1 to 3 weeks, sharing the characteristics of the bursts of LFE (Fig. 5, A and B). The moment-duration relationship of the simulated events occupies a wide space delineated by linear and cubic scaling (Fig. 5C), compatible with various assessments of moment-duration scaling of natural slow-slip events in other tectonic settings (*27*, *28*). However, the moment-duration scaling differs for slow and fast events, implying different source properties linked with slip velocity and shear heating.

The concordance of the model with various seismo-geodetic constraints on the recurrence of LFE swarms and on the underlying slow-slip events gives us confidence to discuss other predictions of the model that cannot be directly assessed based on geophysical data. The evolution of temperature in the unstable temperature-weakening region (Fig. 4) reveals that shear heating during slow-slip events may increase the temperature on the fault by hundreds of degrees, systematically exceeding the solidus of wet granite, i.e., about 600°C at 20-km depth, yet rarely reaching the wet liquidus, around 900° to 1100°C. The distribution of temperature is highly heterogeneous, with the highest temperatures only attained at the center of the slip distributions, where both slip and slip velocities reach their local maximum. During these week-long pulses of high temperature, the maximum velocity attains only about 10^{−7} m/s, far below the seismic regime.

## DISCUSSION

Our results indicate that partial melting may occur around faults due to shear heating, even though the deformation is mostly aseismic, characterized by peak slip velocities only hundred times higher than the tectonic rate (about 1 nm/s) yet much lower than in the seismic regime (about 1 m/s). These inferences are compatible with observations of pseudotachylyte vein injections (*29*) and cockade structures (*30*, *31*) at exhumed fault zones that associate the presence of hot fluids or melts with rapid, localized deformation. However, our results suggest that these and other high-temperature proxies (*29*) may be compatible with the slow-slip regime, reconciling similar observations from laboratory experiments on granitoid cataclasites (*32*).

Our model does not explicitly resolve the seismic emissions, which operate at smaller length scales and shorter time scales. However, given the near-liquidus temperatures reached at some locations, fast slip may be caused by flash weakening associated with partial melting at isolated hotspots. Laboratory experiments at high slip speed indicate that strong weakening may be associated with lubrication of the fault surface by partial melting (*33*), by gouge amorphization (*34*), or by pressurization of pore fluids (*35*), all thermally activated phenomena associated with shear heating. Considering the probable presence of pseudotachylytes, we speculate that flash weakening may be activated during the Parkfield slow-slip events at the locations where near-liquidus temperature is reached, despite the overall low slip velocity. These partial melt pockets may populate the fault at various locations during the seismic cycle following the stress shadows from previous ruptures or be associated with heterogeneity enhancing local shear heating, such as increased normal stress or a thinner active shear zone. Because the regions reaching the highest temperatures represent a small fraction of the slip distribution of any slow-slip event, this may explain the dominance of aseismic deformation in the cumulative moment of slow slip.

Spontaneously occurring deformation in the mid-crust potentially associated with partial melting argues against a top-down control on crustal dynamics (*13*). Instead, the positive feedback between shear heating and temperature-weakening friction provides a mechanism for the spontaneous development of thermal instabilities in the middle crust, even outside any perturbation from the seismic cycle in the seismogenic zone. The complexity of slow-slip event recurrence patterns associated with full and partial ruptures of the temperature-weakening region exemplifies how transient slow slip may redistribute stress and trigger LFE at various times and places, providing an explanatory context for the along-strike variation of LFE activity along the San Andreas fault. If creep is accelerated by weakening processes, the creep-mediated stress transfer may operate efficiently at larger distances than with static or wave-mediated stress transfer.

The collective seismological and geodetic observations at the Parkfield segment of the San Andreas fault illuminate a cohesive picture of the fault zone behavior in the brittle layer (Fig. 6). Considering the velocity dependence of frictional resistance helps define the extent of the seismogenic zone, presumably associated with a narrow cataclasite interface (*16*). Below the seismogenic zone, the fault fabric may consist in a thicker active shear layer made of a distributed network of primary slip surfaces within a fault core. In this active shear zone, the frictional resistance decreases with more elevated temperature and frictional instabilities may develop, depending on the internal configuration of the shear zone, which affects the stability condition. The shear zone may become mylonitic at greater depths, below the brittle-ductile transition.

Thermal instabilities represent a robust interpretation for the slow-slip events in the continental crust, below the seismogenic zone, because the temperature-weakening properties of quartzofeldspathic rocks at these conditions are well constrained from laboratory experiments. Thermal instabilities may occur in other tectonic settings as simultaneous velocity-strengthening, temperature-weakening friction has been inferred for carbonate (*36*), phyllosilicate (*37*, *38*), calcite (*39*), and serpentinite gouges in various hydrothermal conditions. However, slow-slip instabilities below the seismogenic zone at subduction zones may be caused by different processes due to velocity-weakening, temperature-strengthening friction of serpentinite (*40*) in these conditions and the predominance of fluid circulation.

While the velocity-weakening and temperature-strengthening properties of granitic rocks at low temperature, hydrous conditions are important to define the boundaries of the seismogenic zone, the velocity-strengthening and temperature-weakening properties of quartzofeldspathic rocks in shear zones at more elevated temperatures may play an important role on the strength of the mid-crust in the brittle regime and for the rapid redistribution of stress through slow-slip events and slow earthquakes below the seismogenic zone.

## MATERIALS AND METHODS

### Constitutive model for fault slip with thermal activation

We adopt a microphysical model where the slip velocity is controlled by the area of the contact junctions that support the shear and normal loads (*20*), based on the constitutive relationship

(4)where *V* is the sliding velocity, τ is the norm of the shear traction resolved on the fault plane, *T* is the local absolute temperature, μ_{0} is the static coefficient of friction at velocity *V*_{0} and temperature *T*_{0}, A is the real area of contact divided by the nominal surface area, χ is the indentation hardness of the fault material, and *a* ≪ 1 is a power exponent. Fault slip is thermally activated following an Arrhenius law with the activation energy *Q* and the universal gas constant *R*.

The real area of contact depends primarily on the effective normal stress and is modulated by changes in the shape of the microasperities and the quality of contact. Flattening of the microasperities increases the real area of contact junctions and increases the strength of the fault. Comminution during fault slip reduces the average contact area and weakens the fault. Other processes can increase the quality of contact, which may be captured considering an effective contact area. These effects may be captured by a dependence of the real area of contact on effective normal stress and a state variable

$$A=\frac{{\mathrm{\mu}}_{0}\overline{\mathrm{\sigma}}}{\mathrm{\chi}}{\left(\frac{\mathrm{\theta}{V}_{0}}{L}\right)}^{\frac{b}{{\mathrm{\mu}}_{0}}}$$(5)where

$\overline{\mathrm{\sigma}}$ is the effective normal stress accounting for the effect of pore pressure, *L* is the characteristic weakening distance, *b* ≪ 1 is a power exponent, and θ represents the age of contact (*20*).

The constitutive relationship used in the study comes about by first combining Eqs. 4 and 5 to obtain

$$V={V}_{0}{\left(\frac{\mathrm{\tau}}{{\mathrm{\mu}}_{0}\overline{\mathrm{\sigma}}}\right)}^{\frac{{\mathrm{\mu}}_{0}}{a}}{\left(\frac{\mathrm{\theta}{\mathrm{V}}_{0}}{L}\right)}^{\frac{-b}{a}}\text{exp}[-\frac{Q}{R}(\frac{1}{T}-\frac{1}{{T}_{0}}\left)\right]$$(6)then by recasting Eq. 6 to express the frictional resistance as a function of the other controlling variables, to obtain the reciprocal relationship

$$\mathrm{\tau}={\mathrm{\mu}}_{0}\overline{\mathrm{\sigma}}{\left(\frac{V}{{V}_{0}}\right)}^{\frac{a}{{\mathrm{\mu}}_{0}}}{\left(\frac{\mathrm{\theta}{\mathrm{V}}_{0}}{L}\right)}^{\frac{b}{{\mathrm{\mu}}_{0}}}\text{exp}\left[\frac{a}{{\mathrm{\mu}}_{0}}\frac{Q}{R}\right(\frac{1}{T}-\frac{1}{{T}_{0}}\left)\right]$$(7)which is Eq. 1. The parameters *a*, *b*, and *L* serve the same function as in the classical formulation of rate- and state-dependent friction (*21*, *15*), which constitutes the linear terms of a Taylor series expansion of Eq. 7.

Fault strength is modulated by the evolution of the real area and quality of contact through the proposed relationship (*20*)

(8)where *H* is the activation enthalpy for contact healing (*15*). At steady state, the frictional resistance reduces to

(9)indicating that velocity-strengthening behavior is obtained for *a* − *b* > 0, and temperature weakening for *aQ* − *bH* > 0. At steady-state, at the reference velocity *V*_{0} and reference temperature *T*_{0}, the frictional resistance simplifies to

, i.e., the simplified relationship often used as a plastic yield criterion (*20*).

### Dynamics of a spring-slider assembly

We consider the dynamics of a spring-slider assembly with the following governing equation for stress evolution

$$\stackrel{\u0307}{\mathrm{\tau}}=-k(V-{V}_{\mathrm{L}})-\frac{G}{2{V}_{\mathrm{s}}}\stackrel{\u0307}{V}$$(10)where *k* is the elastic stiffness of the spring element, *V*_{L} is the loading rate, *G* is the rigidity, *V* is the sliding velocity, and V_{S} is the shear wave speed. The second term on the right-hand side of Eq. 10 corresponds to the radiation-damping approximation. Eqs. 3, 7, 8, and 10 form a complete set of governing equations that can be solved with the Runge-Kutta method.

A spring-slider assembly is a simple approximation for the dynamics of an isolated fault asperity. We consider the case of a circular velocity-strengthening, temperature-weakening asperity of radius *R* surrounded by rocks of rigidity *G*. Following Eshelby’s stress drop for a circular asperity, the elastic stiffness of the spring element is given by

(11)Assuming *R* = 500 m, *G* = 32 GPa, we obtain *k* = 44 MPa/m. We use *a* = 10^{−2} and *a* − *b* = 4 × 10^{−3}, compatible with the postseismic geodetic measurements following the 2004 Mw 6.0 event (*41*). We drive the slider with a background rate of *V*_{L} = 31 mm/year, i.e., *V*_{L} = 1 nm/s, compatible with the long-term deep creeping rate of the San Andreas fault (*42*). For the reference temperature, we use the steady-state temperature of the evolution law (Eq. 3)

(12)Laboratory experiments provide constraints on the activation energy *Q* and activation enthalpy *H* in the range of 40 to 170 kJ/mol (*15*, *43*), with an average value of 89 ± 23 kJ/mol for wet quartz gouge (*15*). The model shown in Fig. 2 uses *Q* = 90 kJ/mol and *H* = 60 kJ/mol, corresponding to *aQ* − *bH* > 0, temperature-weakening behavior.

Thermal evolution (Eq. 3) depends of the thickness of shear zone *w*, the thermal diffusivity *D*, and the thickness of damage zone *W*. We use *D* = 0.7 mm^{2}/s compatible with laboratory measurements for mid-crustal conditions (*44*). We identify our preferred values for *w*, *W*, and *L* by a Monte Carlo search that penalizes simulated slow-slip events with recurrence times and event duration significantly different from 1 year and 1 month, respectively. For each set of thermal-mechanical model parameters, we simulate 300 years of fault dynamics. After discarding the first 200 years of simulation to mitigate the undesirable effects of initial conditions, we estimate the duration and recurrence times of slow-slip event. We explore the model space with about 30,000 such forward models. The Monte Carlo sampling provides an estimate of the joint probability density function of the model parameters. Figure S1 shows the marginal distribution of the active shear zone thickness *w*, the diffusion distance *W*, and the characteristic weakening distance *L*. Figure S1 also showcases the marginal bivariate probability densities to represent the tradeoffs between model parameters. Reasonable models that produce slow-slip events with acceptable recurrence times and duration can be obtained for a realistic range of model parameters (fig. S2), with 10^{−4} < *w* < 10^{−3} m, 1 < *W* < 10 m, and 1 < *L* < 10 mm. The model presented in Fig. 2 uses *W* = 1.7 m, *w* = 0.14 mm, and *L* = 3 mm. This combination of physical properties reconciles laboratory constraints and provides a reasonable range of temperature, peak velocity, and recurrence time representative of the shallow LFE families northwest of Parkfield and their associated slow-slip events. The other relevant variables are provided in table S1.

### Modeling thermal instabilities on a planar fault

We consider a planar fault embedded in an elastic half-space with the following governing equation for stress evolution

$$\stackrel{\u0307}{\mathrm{\tau}}(\mathbf{x})={\int}_{\mathrm{\partial}\mathrm{\Omega}}\mathbf{K}(\mathbf{x};\mathbf{y})\xb7(\mathbf{V}(\mathbf{y})-{\mathbf{V}}_{\mathrm{L}})d\mathbf{y}-\frac{G}{2{V}_{\mathrm{S}}}\stackrel{\u0307}{\mathbf{V}}(\mathbf{x})$$(13)where **K**(**x**;**y**) is stress kernel relating the slip velocity at position **y** to the traction on the fault at position **x**, **V** is the instantaneous sliding velocity vector encompassing the strike-slip and dip-slip components, and *V*_{L} is the loading rate vector with a nontrivial component only in the strike direction. The fault area is denoted by ∂Ω. The numerical simulation is accomplished using the Unicycle software (*40*, *8*), based on the fourth/fifth-order Runge-Kutta solver and the boundary-integral method.

We resolve fault dynamics on a fault segment extending 45 km in length and 16 km in width, starting at 10 km depth, 12 km northwest of the epicenter of the 2004 Mw 6.0 Parkfield earthquake (Figs. 1C and 3B). The fault is vertical, striking 39° to the northwest. The fault plane is discretized with 200 m square patches, sufficient to capture the slow nucleation of thermal instabilities. All the simulated domain is velocity strengthening and temperature weakening. However, we include a 15-km-long, 1-km-wide, unstable region centered on the shallow LFE family with a thin shear zone thickness of *w* = 0.14 mm, promoting efficient shear heating. The region surrounding this patch is stable, with *w* ≫ 1 mm, which inhibits efficient shear heating. The model shown in Figs. 2 and 3 uses *w* = 1 m for the stable region, but we found that any value larger than 1 mm would produce similar results. The unstable region is also centered in the geodetically inferred location for the correlative slow-slip events (*11*).

For the background temperature, we investigate different profiles for the continental crust based on the thermal regime of the Parkfield crust (*45*). First, we assume a uniform background temperature of *T*_{b} = 400°C, a value representative of the temperature found at 17-km depth, in the narrow 1-km-wide unstable region. Second, we consider a thermal gradient of 21°C/km, leading to a temperature of 372°C in the middle of the temperature-weakening zone. The temperature at depths shallower than 10 km and deeper than 26 km, i.e., outside the modeled domain, is not relevant. We assume an effective normal stress of 20 MPa and a shear modulus of *G* = 40 GPa (*46*). We apply the long-term loading rate *V*_{L} = 10^{−9} m/s, compatible with the geodetically determined tectonic loading rate of 31 mm/year. The other parameters are similar to those adopted in the spring-slider model and are summarized in table S2. Representative simulation results for the model with a geothermal gradient are shown in Figs. 3 and 4. A comparison of the statistics of recurrence time and duration of slow-slip events for the uniform-temperature and thermal-gradient models is shown in fig. S3. The choice of the background temperature in the stable temperature-strengthening region has little impact on the model outcome because of the narrow depth range of the unstable region.

### Duration and recurrence time of LFE bursts

The seismic activity that is clearly associated with transient slow-slip includes a cluster of 10 LFE families (yellow circles with black contour in Fig. 1). For each of the 10 families, we group the densely clustered tremors into different bursts. The sequential bursts are separated by a quiescent period of days. The resulting LFE bursts plotted against time is shown in Fig. 1B. We select the bursts that include more than 60 LFE occurrences following the same criterion used for geodetically detecting slow-slip event (*11*). We quantify the quiescent time between the so-defined bursts. Their recurrence frequency and a comparison with the simulated slow-slip events are shown in Fig. 4E. A comparison of the simulations assuming a uniform background temperature or a finite thermal gradient is shown in fig. S3.

**Acknowledgments: **We appreciate the comments from three anonymous reviewers. **Funding:** L.W. was supported by the National Natural Science Foundation of China under grant numbers NSFC-41674067 and NSFC-U1839211. S.B. was supported in part by the National Science Foundation under award number EAR-1848192. **Author contributions:** S.B. designed the study. L.W. and S.B. conducted the study and wrote the manuscript. **Competing interests:** The authors declare that they have no competing interests. **Data and materials availability:** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).