Research ArticleGEOPHYSICS

Excitation of San Andreas tremors by thermal instabilities below the seismogenic zone

See allHide authors and affiliations

Science Advances  04 Sep 2020:
Vol. 6, no. 36, eabb2057
DOI: 10.1126/sciadv.abb2057


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.


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.

Fig. 1 Deep LFEs and slow slip along the San Andreas fault.

(A) LFE (10) (yellow circles) and Mw ≥ 3 seismicity from the Northern California Earthquake Data Center (NCEDC) catalog (gray circles) for the period 2004–2017. The 10 LFE families above 20 km depth are marked with a black contour. The Mw 6.0 Parkfield earthquakes of 1966 and 2004 are marked by the white and red stars, respectively. The displacements (black vectors) at GPS sites (white triangles) are due to a simulated large slow-slip event of Mw 4.72, whose slip distribution is shown in Figure 3B. (B) Temporal behavior of the shallowest 10 LFE families northwest of Parkfield. (C) Vertical cross section along the San Andreas fault. The background color, associated with a correlation coefficient, indicates the likely location of a slow-slip event driving the shallowest LFE (12).

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.


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τ=μ0σ¯(VV0)aμ0(θV0L)bμ0exp[aμ0QR(1T1T0)](1)where τ and σ¯ are the shear and effective normal stress on the fault; V and θ are slip velocity and age of asperity contact, respectively; μ0, V0, and T0 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θ̇=exp[HR(1T1T0)]VθL(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 = T0, 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Ṫ=DW2(TTb)+τVwρc(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, Tb 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).

Fig. 2 Thermal instabilities in a spring-slider assembly.

(A) Schematic strength profile showing the velocity-weakening shallow crust above a velocity-strengthening mid-crust and a ductile lower crust. A velocity-strengthening, temperature-weakening fault patch undergoes shear heating in the active shear layer of thickness w. Temperature diffuses to the bath temperature Tb through a fault zone of thickness W. (B) Dependence of steady-state friction on velocity and temperature. (C) Repeat cycles among velocity and temperature, and shear stress and temperature. (D) Repeat cycles among age of contact and temperature, and velocity and temperature. Different phases are indexed ordinally and labeled with different colors. (E) Evolution of velocity and temperature during three slow-slip events. (F) Corresponding evolution of the age of contact and shear stress.

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 ab < 0 and temperature strengthening for aQbH < 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 (ab > 0) but temperature-weakening (aQbH > 0) fault patch embedded in a rock matrix at temperature Tb = 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.

Fig. 3 Numerical simulations of slow-slip events on the San Andreas fault.

(A) Surface displacements produced by two slow-slip events representing partial ruptures of the unstable temperature-weakening area. (B) Slip distribution of slow-slip events with contour lines for cumulative slip of 10, 20, and 30 mm. The black contour is for the slow-slip event of Mw 4.72 shown in Fig. 1. (C) Simulated fault-parallel surface displacements for GPS stations P292, P293, P289, and P291. (D) Slip history at the center of the unstable patch. (D) Temporal evolution of slip velocity (black) and temperature (gray) at the patch center.

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).

Fig. 4 Spatiotemporal evolution of slip velocity and temperature along the San Andreas fault.

(A) Horizontal velocity profile A-A′ at 17-km depth showing a complex sequence of slow-slip events. The inset locates the profile. Light gray and dark lines indicate contour of velocity of 1 and 3 nm/s, respectively. Dark profiles indicate velocity above 3 nm/s. (B) Evolution of maximum velocity anywhere on the fault. (C) Evolution of temperature on profile A-A′. The contours are the same as in (A). (D) Evolution of maximum temperature anywhere on the fault.

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.

Fig. 5 Recurrence pattern of Parkfield LFE and simulated events.

(A) Histograms of recurrence times for simulated slow-slip events (bars) and for the 10 shallow LFE sources (thick profile). (B) Histograms of event duration for simulated slow-slip events and the same LFE sources. (C) Moment-duration scaling for the catalog of simulated slow-slip events. The scaling for linear (solid line) and cubic (dashed line) relationships is shown for reference. Two clusters of events can be determined based on the peak velocity, below (squares) or above (triangles) Vpeak = 1.6 × 10−8 m/s, associated with different moment/duration scaling.

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.


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.

Fig. 6 Schematic of mechanical properties of the San Andreas fault in the brittle regime.

(A) Depth dependence of the steady-state frictional properties ab that control the velocity dependence. (B) Depth dependence of the steady-state frictional properties aQbH that control the temperature dependence. Thermal instabilities may occur in a temperature-weakening domain, depending on stability conditions. (C) Fault fabric at depth and seismic cycle behavior. (D) Horizontal section of a shear zone at the root of the San Andreas fault with multistranded slip surfaces. (E) Details of the fault core with primary slip surfaces associated with shear heating and pseudotachylyte injection veins resulting from partial melting during slow-slip events.

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.


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 relationshipV=V0(τAχ)μ0aexp[QR(1T1T0)](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 V0 and temperature T0, 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 variableA=μ0σ¯χ(θV0L)bμ0(5)where σ¯ 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 obtainV=V0(τμ0σ¯)μ0a(θV0L)baexp[QR(1T1T0)](6)then by recasting Eq. 6 to express the frictional resistance as a function of the other controlling variables, to obtain the reciprocal relationshipτ=μ0σ¯(VV0)aμ0(θV0L)bμ0exp[aμ0QR(1T1T0)](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)θ̇=exp[HR(1T1T0)]VθL(8)where H is the activation enthalpy for contact healing (15). At steady state, the frictional resistance reduces toτ=μ0σ¯(VV0)abμ0exp[aQbHμ0R(1T1T0)](9)indicating that velocity-strengthening behavior is obtained for ab > 0, and temperature weakening for aQbH > 0. At steady-state, at the reference velocity V0 and reference temperature T0, the frictional resistance simplifies to τ=μ0σ¯, 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τ̇=k(VVL)G2VsV̇(10)where k is the elastic stiffness of the spring element, VL is the loading rate, G is the rigidity, V is the sliding velocity, and VS 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 byk=7π16GR(11)Assuming R = 500 m, G = 32 GPa, we obtain k = 44 MPa/m. We use a = 10−2 and ab = 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 VL = 31 mm/year, i.e., VL = 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)T0=Tb+μ0σ¯VLwρcW2D(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 aQbH > 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 mm2/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τ̇(x)=ΩK(x;y)·(V(y)VL)dyG2VSV̇(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 VL 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 Tb = 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 VL = 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.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


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.

Stay Connected to Science Advances

Navigate This Article