Dynamic weakening during earthquakes controlled by fluid thermodynamics

Earthquakes result from weakening of faults (transient decrease in friction) during co-seismic slip. Dry faults weaken due to degradation of fault asperities by frictional heating (e.g. flash heating). In the presence of fluids, theoretical models predict faults to weaken by thermal pressurization of fault fluid. However, experimental evidence of rock/fluid interactions during dynamic rupture under realistic stress conditions remains poorly documented. Here we demonstrate that the relative contribution of thermal pressurization and flash heating to fault weakening depends on fluid thermodynamic properties. Our dynamic records of laboratory earthquakes demonstrate that flash heating drives strength loss under dry and low (1 MPa) fluid pressure conditions. Conversely, flash heating is inhibited at high fluid pressure (25 MPa) because water’s liquid–supercritical phase transition buffers frictional heat. Our results are supported by flash-heating theory modified for pressurized fluids and by numerical modelling of thermal pressurization. The heat buffer effect has maximum efficiency at mid-crustal depths (~2–5 km), where many anthropogenic earthquakes nucleate.

Theoretical models predict that TP dominates at large slips and/or above mid-crustal depths, while FH is the dominant weakening mechanism at small slips and/or greater depths 12 . However, these models rarely incorporate water thermodynamics (notably phase transitions). Moreover, in contrast to FH, TP during seismic slip has not been established experimentally as of yet 8,17 , due to the difficulty of reproducing spontaneous dynamic rupture under realistic stress and fluid pressure conditions in the laboratory. In summary, fully controlled experiments studying dynamic shear instabilities in the presence of fluid pressures (pf) have been lacking so far and hydro-thermo-mechanical couplings during dynamic rupture remain still unclear.
Here we conducted stick-slip experiments (laboratory proxies for earthquakes 19 ) on Westerly Granite (WG) saw-cut samples (Supplementary Fig. 1) under triaxial stress conditions (principal stresses σ 1 > σ 2 = σ 3 ). The experiments were done at stresses representative of the upper continental crust 19 (effective confining pressures σ 3 ' = σ 3 − pf = 70 MPa). We imposed three different fluid pressure levels: dry, low fluid pressure and high fluid pressure, which correspond to pf = 0, 1, and 25 MPa, respectively (hereafter referred to as dry, Low Pf and High Pf , respectively). Combining dynamic stress evolutions with on-fault resolved displacements and microstructural analysis of the postmortem specimens evidenced that distinct dynamic weakening mechanisms (FH and TP) were activated at the different fluid pressure levels. Further, we applied thermal weakening models to our experimental data including the evolution of fluids thermophysical properties with pressure and temperature. The results showed that both FH and TP were activated during co-seismic slip and that their relative contributions are controlled by the evolution of water's thermophysical properties at phase transitions.

Results and discussion
Mechanical results. Continuous records of shear stress (τ) versus time (Fig. 1a) showed that, during a stick-slip event, τ dropped from an initial peak static value (τ 0 ) down to a final residual value (τ f ), resulting in a static stress drop (Δτ s = τ 0 − τ f ). Highfrequency records 9,10 ( Fig. 1b-d) showed that τ first increased from τ 0 up to a peak dynamic value τ p and then abruptly dropped to a minimum dynamic value τ min before recovering to τ f , thereby defining a dynamic (or breakdown) stress drop (Δτ b = τ p − τ min ). The dynamic rise of τ up to τ p resulted from stress amplification (i.e. stress intensity factor) at the rupture tip 20 .
Earthquakes during dry experiments presented larger Δτ s values ( Fig. 2a) (from 30 to 45 MPa versus from 10 to 30 MPa and from 5 to 18 MPa under Low Pf and High Pf , respectively), while larger breakdown stress drops were recorded under Low Pf . There Δτ b was 14% larger on average than in dry conditions and was remarkably 73% higher than at High Pf . Note that τ 0 (i.e. the amount of elastic energy stored in the system) was similar at both fluid pressures but smaller than in dry conditions. The total slip per event was similar for dry and Low Pf conditions but was two times smaller under High Pf . In all conditions, peak static frictional strength (μ 0 = τ 0 /σ n0 ') ranged between 0.6 and 0.89 (Fig. 2b, compatible with Byerlee's law 21 ) but was approximately 17% higher in dry experiments than at Low Pf and High Pf . This indicates lower static shear strengths in the presence of fluids that resulted from a reduction of adhesion along fault surface in the presence of water 22 (i.e. to a decrease of the contact surface energy). Regarding weakening processes, the dynamic friction (μ d = τ min /σ n ', see Methods) was lower at Low Pf (from 0.02 to 0.24) than at dry (from 0.29 to 0.39) conditions and High Pf (from 0.42 to 0.51). Such differences in the dynamic fault strength (i.e. evolution of τ and μ d ) in these three experiments (performed at constant σ 3 ') suggest the activation of distinct weakening mechanisms during earthquake rupture. Such mechanisms seem less effective at High Pf (i.e. smaller slip and higher μ d for equivalent τ 0 and Δτ s ) and slightly more effective at Low Pf (i.e. larger Δτ b and lower μ d , leading to a transient, almost total strength loss).
Microstructural observations. Scanning electron microscopy on postmortem fault surfaces (Fig. 3, Supplementary Fig. 4) revealed 20 μm long patches of ropey, stretched material elongated along the shear sense in both the dry (Fig. 3b) and Low Pf (Fig. 3c) experiments. These structures are consistent with melting of fault asperities during co-seismic slip and may explain the strong weakening observed in these experiments 9,10,18 . Conversely, such structures were not found at High Pf (Fig. 3d) where the surfaces  Each curve corresponds to one stick-slip event, change in colour hue accounts for different events. In addition to τ 0 and τ f , the maximum and minimum dynamic values of shear stress, τ p and τ min , are presented as examples, defining a breakdown stress drop (Δτ b = τ p − τ min ). b Dry experiment, red curves. c Low fluid pressure experiment, blue curves. d High fluid pressure experiment, black/grey curves were covered with asperity debris of sizes~0.2-5 μm. No evidence of melt was found, suggesting that lower temperatures were reached at the asperity contacts and confirming the reduced efficiency of frictional heating and weakening at High pf .
Flash temperature in the presence of pressurized fluids. To support our experimental and microstructural results, we computed the flash temperatures (maximal transient temperatures) reached on 20 μm radius asperities due to shear heating 2,23 during their lifetimes (t c ) (Methods, Fig. 4a). In the presence of fluids, water cools asperities through heat capacity and latent heat (acting as a heat barrier) of a finite water volume surrounding the highly stressed asperity 7 (Fig. 4b). The main hypothesis of the model is that the fluid volume surrounding the asperities is at thermal equilibrium with the asperity. This assumption should remain valid during frictional slip since the thermal diffusion length ( ffiffiffiffiffiffiffiffiffiffiffiffi ffi π:κ th :t p with κ th the thermal diffusivity and t the heating time) is close to the asperity size at FH velocities (commonly admitted~10 cm s −1 (refs. 2-7 )), and when the solid-solid contact starts slipping, a liquid-solid contact forms immediately, allowing for fast temperature equilibrium between the asperity and the surrounding fluid. Conversely to previous studies 7, 24 , we also included the isobaric evolution of water's specific and latent heat (c pw and L w ), as well as density (ρ w ) with temperature in the calculation 25 (Fig. 4c-e, Methods, Supplementary Fig.7). Given our experimental conditions, we considered water cooling of asperities as a purely diffusive mechanism (no advection) for fault permeabilities <10 -17 m 2 at Low pf and <10 -18 m 2 at High pf (Methods, Supplementary Fig. 5). Note that in this model we considered the maximum temperature that can be reached at asperities affected by the cooling effect of water 7 . Such temperature differs from the temperature history at asperities during slip. Under dry conditions, when no buffering takes place, the flash temperature becomes the exact solution for the onedimensional heat diffusion problem 2,10,26 at the asperity scale (assuming the contact shear stress at asperities rather than the macroscopic shear stress distributed along the interface). There temperature increased as a power law of slip (see Supplementary  Fig. 6 for other asperity sizes). The FH temperature (approximately 1000°C 2-5 ) was reached for slip rates >10 cm s −1 during the asperity lifetime, as predicted by FH theories and previous experiments [2][3][4][5][6][7] . At those velocities, in the Low Pf case, waterbuffered temperatures were observed in the first half of the contact lifetime, and so, flash temperatures remained <179°C, i.e. while water stayed in a liquid state. Longer slips at such seismic velocities (e.g. higher frictional heat) allowed water to overcome the liquid-vapour phase transition temperature during t c, inducing a strong drop in ρ w and c pw (roughly falling to 0.5% and 50% of their room temperature values, respectively; Fig. 4c, d). In this case, vaporization of water enhanced shear heating at contacts and allowed FH of asperities for slip velocities larger than~10 cm s −1 , as also observed in dry conditions. Conversely, at fluid pressures ranging from 25 to 70 MPa, temperature rise was strongly buffered by water cooling during t c due to the liquid-supercritical transition. This phase change requires a distributed amount of energy over a finite temperature range, opposed to the case of isothermal vaporization where L w acts as a heat barrier. Therefore, the heat capacity of water increases by 1400% during the transition at pf = 25 MPa (Fig. 4d) while the drop in density is smoother than in the case of vaporization 25 (Fig. 4c). At high fluid pressures, water turned out to be an extremely efficient heat buffer, inhibiting FH phenomena and hindering rises in temperature to the liquid-supercritical phase transition temperature (~373°C at pf = 25 MPa, Supplementary  Fig. 7) at asperity contacts during their lifetime. Temperature rise was buffered even for slip rates of 1 m s −1 (admitted slip rate during earthquake propagation 1,2,4,6 ). This major heat sink explains the reduced dynamic weakening observed at High Pf and the absence of frictional melt on the fault surfaces.
Shear heating and TP of fluid saturated faults. The liquid-vapour transition has been thought to have strong thermal effects on faulting, inhibiting temperature rise due to TP during co-seismic slip 14,15 . In high-velocity friction experiments 8,14 , TP enhanced the friction drop of~0.1, which is comparable to the difference observed between the dynamic friction recorded during Low Pf and dry conditions (Fig. 2b). Such difference could also be due to a reduction of melt viscosity through hydration in the presence of fluids. However, rotary shear experiments have demonstrated that the chemical compositions of melts developed after long slip times (>10 s) under vacuum, room humidity and fluid-saturated conditions were identical 7 , discarding the possibility of melt-hydration in our experiments (here the total slip time was <30 μs). TP due to fluid pressurization could then be a candidate to explain the slightly lower dynamic friction values observed at Low pf while FH remains the dominant weakening mechanism.
While FH explains the large dynamic strength drops observed in dry and Low Pf conditions, it does not explain the small strength drops observed at High Pf conditions. To explain the . b Schematic top view of the considered contact geometry, in black we see the asperity contact of radius r = 20 μm stressed at a normal stress σ c and at a shear stress τ c . The asperity is surrounded by a volume of water V w , which buffers temperature. c Temperature versus water density 25 for the different fluid pressures in MPa. d Temperature versus water specific heat 25 . e Temperature versus water's latent heat 25 . Note that the latent heat (L w ) decreases with rising pressure until it vanishes at the critical point (~22 MPa and~373°C) small stress drops observed at High Pf , we computed the temperature evolution on a bulk planar fault in both drained (Fig. 5a) and undrained (Fig. 5b, c) conditions using a finite difference numerical model (e.g. Methods). We considered full thermodynamic evolution of fluid properties with pressure and temperature 15,16,25 . Under drained conditions, we observed that the reached temperatures (which are a maximum estimation of the possible temperature in the bulk fault since the shear stress for heat generation is taken as the static fault's shear strength of our experiments) remained below rock's thermal degradation temperature (~1000°C) even for slips larger than the maximum slip observed in the experiments (~250 μm) (Fig. 5a). This observation is in agreement with our microstructural observations, since melting was not pervasive over the sample surface but was localized at asperity scale ( Fig. 3c and Supplementary Fig. 4c), as predicted by FH theory. Nevertheless, under our experimental conditions, the water heat buffer effect due to the liquid-supercritical transition was still observed for initial pore pressures >22 MPa. Note that strong temperature rises on the fault due to instantaneous water vaporization took place in a similar manner than for the flash temperature computations, confirming our calculations at the asperity scale. Under undrained conditions (Fig. 5b, c), since the fault's stress obeys the effective pressure law, we observed an initial fast decay in friction due to TP. The decay then stabilized leading to friction drops of~0.1-0.15 for slips of~20-150 μm in all fluid pressure conditions. Such friction drop values are remarkably consistent with the friction drops observed in High Pf experiments (Fig. 2b). Therefore, at High Pf , TP might well be the dominant weakening mechanism in our experiments, as supported by our microstructural analysis.
Implications for natural and induced earthquakes. Similar stress evolutions observed in experiments conducted at other effective stresses and at pf = 45 MPa ( Supplementary Fig. 3) suggest that the observed heat buffer operates even at higher fluid pressures, where the liquid-supercritical transition is smoother 25 (e.g. Fig. 4c, d and Supplementary Fig. 7). To further study the depth dependence of this heat buffer effect, we computed again the temperature rises (in both drained (Fig. 6a) and undrained conditions (Fig. 6b, c)) due to TP with a depth extrapolation. The extrapolation was done for a mean stress equal to the lithostatic overburden gradient of 27 MPa km −1 , a hydrostatic fluid pressure rise of 10 MPa km −1 , a geothermal gradient of 30°C km −1 (ref. 27 ) and an initial friction of 0.7 (e.g. Methods). We observed that a heat buffer can operate for fluid pressures up to 45 MPa in both drained and undrained cases but its efficiency is strongly reduced when fluid pressures reach 70 MPa (~7 km depth). At large depths, higher background stress and a smoother supercritical transition allow to overcome the transition temperature for smaller slips when sliding at seismic slip rates (~1 m s −1 ), consistent with previous studies on the depth dependence of weakening mechanisms 12 . Nevertheless, the dynamic friction values predicted by TP theory are similar at all depths for a given final slip, likely because at greater depths (>7 km), the background driving stress has a stronger effect than the pore fluid pressure rise on TP at small slip 12 . Previous TP models considered the liquid-supercritical transition and found no significant effect of the transition on dynamic ruptures 16 but did not consider the effect of FH at the microscopic level. Here we demonstrate that water phase transitions may control FH at the asperity level by acting as a major heat buffer and so they can control earthquake rupture. The initial fluid pressure level is a critical parameter that cannot be neglected via the effective pressure concept because it controls water thermodynamics. Thus dependencies on temperature and pressure of thermodynamic fluid properties should be taken into account in future weakening models, in particular at the microscopic level 2,12,13 . Extrapolation of our results to natural conditions suggests that the heat buffer effect has a maximum efficiency at mid-crustal depths (~2-5 km) where major anthropogenic earthquakes appear 28 .  Fig. 1) (40 mm diameter and 88 mm length). This material was selected because it is simple in composition, it is representative of the continental crust and because of its very fine grain (<1 mm), perfect homogeneity, isotropy and low alteration degree. A thermal treatment at atmospheric pressure was performed on the samples in order to increase the samples' crack density (i.e. permeability) and so to allow good saturation and reasonable fluid diffusion times in the samples. Samples were heated at a gradient of 5°C min −1 up to 450°C. Then the target temperature was kept constant for 2 h. Finally, the furnace was turn off and the samples were let to cool down overnight. The temperature ramp imposed in the furnace allowed having no temperature gradient inside the sample during heating. Since the thermal diffusivity of the rock is~κ th = 10 -6 m 2 s −1 , and with a sample radius of r = 2 × 10 -2 m, temperature should equilibrate in r 2 /κ th = 400 s (~7 min). Thermal cracking arose from differential thermal expansion of neighbouring grains thus allowing increasing both intergranular and grain boundary microcracking without overcoming the Quartz Alpha-Beta transition (578°C). Permeability was measured after thermal treatment and was~5 times higher than that of untreated samples (from 2 × 10 -19 m 2 to 1 × 10 -18 m 2 at 5 MPa confining pressure). The cylinders were then cut to the correct length, and the top and bottom bases grinded to ensure perfect planarity with the horizontal. Then the samples were saw-cut at an angle (θ) of 30°to the sample's long axis to create an artificial elliptical fault of major axis L = 80 mm and minor axis l = 40 mm. The apparent contact area being: Fault's surfaces were then grinded to ensure perfect contact in the fault and then roughened with #240 grit paper to ensure a minimum cohesion along the fault's interface and impose a constant fault roughness in all the specimens.
Experimental set-up. The apparatus used was the tri-axial press of ENS Paris built by Sanchez Technologies. It is a servo-controlled oil medium confining cell with maximum confining pressure of 100 MPa. Axial loading was controlled by a separated servo pump acting on an axial piston (maximal stress of 680 MPa on 40 mm diameter samples). Fluid pressure regulation was assured by a double syringe pump (Quizix 20k) able to reach 120 MPa fluid pressures (1 kPa pressure accuracy, 1 μL volume accuracy). Under this configuration, shear stress (τ), normal stress (σ n ) and slip (D f ) resolved on the fault can be expressed as: where σ′ refers to effective stress as σ′ = σ − pf and where ε 1ext is the measured axial strain on the whole system; ε 1s is the axial strain of the sample corrected by the stiffness of the apparatus using linear elasticity, Δσ = (σ 1 − σ 3 ) the differential stress; E ap the stiffness of the apparatus; L the sample's length; D 1 the axial displacement of the sample and, finally, D f the projected displacement on the fault. Finally, making the reasonable assumption that confining pressure (σ 3 ) does not change during stick-slip events, near fault friction (μ) is calculated as: The recorded parameters during deformation were as follows. In the far-field, we recorded the axial and confining pressures through pressure transducers of 0.001 MPa resolution. In addition, axial displacement was measured by recording three Foucault current sensors of 0.1 μm resolution. The sampling rate on far field sensors was 100 Hz. These provided the macroscopic deformation of the system (sample plus apparatus deformation). In the near-field, we measured stress and strain through strain gages glued 3 mm away from the fault (Supplementary Fig. 1). Gages were coated with a cyanoacrylate gel, which prevented shortages due to pressurized water. These sensors allowed a local recording of the principal strains at 10 MHz sampling frequency. A full (Wheatstone) bridge configuration gage (HBM 3/350 VY41) allowed measuring directly the differential strain (ε 1 − ε 3 ). To calibrate the gage, we measured the constant Young modulus of the rock during elastic loading phase. Then we had direct conversion from the strain recorded at the gage to the corresponding far field measured stress using the measured sample's Young modulus (differential stress (σ 1 − σ 3 )) ( Supplementary Fig. 2). Gages allowed to record the dynamic stress change of each stick-slip event through an acoustic emission trigger set-up 9,10 at 10 MHz sampling frequency. Strain gage data were recorded continuously at 100 Hz to observe the overall evolution of near fault shear stress.
Loading procedure and laboratory earthquakes. Stick slip experiments were performed under nominally dry and fluid pressure conditions. Confining pressures ranged from 50 to 95 MPa, and fluid pressure from 0 (i.e. dry) to 45 MPa. Constant strain rate was imposed at~1 × 10 -5 s −1 (see Supplementary Table 1 for the detailed experimental matrix). The experimental procedure was as follows: We first increased the confining and axial pressures up to 10 MPa. Then, in case of pressurized fluid experiments, we carefully flushed air away from the sample, then we increased the fluid pressure to 5 MPa in both upper and lower reservoirs. We then waited for pressure and volume equilibrium between the two reservoirs. The axial, confining and fluid pressures were increased (fluid pressure was decreased in low fluid pressure experiments) in parallel up to their target values. Again, we waited for fluid pressure and volume equilibrium between the reservoirs. Finally, axial pressure was increased at constant axial loading rate while fluid and confining pressures were held constant. Both shear and normal stresses increased with axial loading until shear stress reached the strength of the fault. At this point, the stickslip instability occurred and was accompanied by a brutal release of shear stress and seismic slip on the fault plane. Such an event is a stick-slip event or laboratory earthquake.
The reproducibility of the shear stress evolution for successive events in each configuration suggests that the possible change in surface topography with increasing number of events did not affect the fundamental processes accounting for earthquake rupture propagation.
Flash temperature computation. Flash temperature is the maximum transient temperature responsible for fast weakening of fault frictional strength during sliding 23 . Flash temperature is reached at an asperity for the lifetime of contacting asperities (t c ) and depends on slip rate v in (m s −1 ), material properties (in particular thermal diffusivity κ th in (m 2 s −1 )) and asperity radius r in (m).
Strong frictional weakening happens when the temperature rise at the contacting asperity reaches values close to the melting or thermal decomposition temperature of the rock, which can be taken equal to 1000°C for many rock lithologies [2][3][4] .
We considered a frictional interface where the real contact area (A r ) is only a fraction of the nominal contact area (A) 29,30 . The major part of the contact is held by asperities that deform mostly plastically and are stressed closely to their yield strength 29,30 . For simplicity, we considered rounded asperities of radius (r) and height (h). In the presence of fluid pressure, we defined a volume of water (V w ) that interacted thermally with the highly stressed asperities such that V w corresponded roughly to the volume of water displaced by the contact sliding during its lifetime. Such volume of water is in convective contact with the asperity and is defined 7 as V w = h.π.((2r) 2 − r 2 ) (Fig. 4b).
To compute the temperature elevation per unit surface at a contacting asperity during slip acceleration, we consider that this elevation is due to diffusion of a heat source rate τ c ×v where τ c is the shear stress at the contact and v is an arbitrarily increasing slip velocity. As argued by Rice 2 , the heat input per unit surface over the contact time (and so, until weakening takes place) is directly related to time of contact defined as: t c = r/v. Then we define the flash temperature rise as a heat input term due to shear at the slip rate v and a temperature buffering term 7 due to the volume of water surrounding the asperities V w as defined in the geometry and shown in Fig. 4b such that: Tflash = f(τ c , ν) − g(T, ρ ω (P, T), c pw (P, T)).
Notice that in this simple model we did not consider the evolution of density or specific heat of solid asperities with temperature (see final equation) and that possible dynamic changes of contact hardness and other material properties due to the flash temperature rise 31 were neglected. Here the fluid pressurized in the fault zone is water. The isobaric evolution of water's specific heat and density with temperature at the experimental pressures were taken from NIST database for thermophysical properties of fluids 24 (based on the IAPWS97 industrial thermodynamic formulation) at different imposed fluid pressures.
The following considerations were used for this model: Asperities of radius r and height h. In a frictional interface, the real contact area of the two surfaces involved is substantially smaller than the apparent contact area 29,30 . Therefore, the load supported by each contacting asperity is considerably higher than the normal stress applied to the apparent surface. In our experiments, microstructural analysis (Fig. 3) showed an initial asperity sizes of~2-40 μm (Fig. 3a). After deformation, the melted patches in dry and low fluid pressure experiments had maximal sizes of~20× 20 μm 2 . Therefore, we defined the maximum asperity size of a radius r = 20 μm and an asperity height h = r = 20 μm.
Applied forces considered: shear and normal stress. Here peak shear stress considered for all simulations matched the average peak static shear stress found during low and high pore pressure experiments, τ 0 = 70 MPa. The peak friction reached was averaged to μ 0 = 0.7 respecting 'Byerlee's rule' 21 . Therefore, the peak static normal stress considered in this model was: σ n0 = τ 0 /μ 0 = 100 MPa. If α = A r / A is the ratio between real contact area and nominal contact area, it writes α = P m / σ n0 29 . Where P m = 6 GPa is the estimated penetration hardness of WG taken as a weighted average 32 of hardnesses of the minerals present in the granite. The shear stress held by a single asperity writes then τ c = α.τ 0 = 4.2 GPa.
Pure diffusion in the vicinity of the contacts surrounded by pressurized fluid. The question arises whether cooling process by convection of water surrounding the asperities is a purely advective, mixed advective/diffusive or purely diffusive process. Owing to the intense and fast heating developed at the highly stressed contacts, the interacting water volume V w was considered instantaneously in thermodynamic equilibrium prior to weakening. The temperature of the fluid is therefore that of the contacts. We then calculated the ratio between the time needed for temperature at asperities to equilibrate with water by advection (t adv ) and the time needed for temperature at asperities to diffuse in water (t heat ) ( Supplementary  Fig. 5). In the same water volume, the ratio: t adv /t heat is expressed as the inverse ratio of hydraulic and thermal diffusivities, respectively. D hy = k(η*.β) −1 is the hydraulic diffusivity where k is the in plane fault's permeability, η* is the fluid viscosity and β is the storage capacity of the interacting volume. D th = λ*(ρ w *.c pw *) −1 is the thermal diffusivity of the fluid volume, where λ* is thermal conductivity of the fluid, ρ w * is the fluid density and c pw * is the fluid specific heat. All values marked with * are dependent on temperature 25 . The higher the t adv /t heat ratio, the more the cooling process is diffusive. On the other hand, for low values of t adv /t heat (<1), the cooling process should be highly enhanced by fluid circulation in the fault and so it becomes an advective process. These calculations showed that the heating process is purely diffusive in the low pressure case (1 MPa) for fault permeabilities <10 -17 m 2 in the whole temperature range. At high fluid pressure (25 MPa), the process is purely diffusive for permeabilities <10 -18 m 2 . At the normal stresses developed in our experiments, in the absence of fault gouge, we estimate that the permeability of the fault was close to that of the surrounding material and so close to values probably inferior to 10 -18 m 2 . We conclude that a purely diffusive model represents well the cooling effect of water during FH in the vicinity of asperities. Thus the model assumed that a finite volume of water in the vicinity of the asperity interacted thermally with the latter through heat capacity and latent heat of vaporization (Fig. 4c-e). Calculations presented in Supplementary Fig. 5 are in agreement with the experimental results while fault permeabilities are >10 -17 m 2 for low fluid pressure. In the case of high fluid pressure experiments, the cooling process should be enhanced by advection around the contacts at the temperature of the liquid/supercritical transition for fault permeabilities reaching 10 -19 m 2 .
Finally, using the stated parameters and considerations, a heat balance per unit area at the asperity where the heat stored in the asperity (V asp .ρ.c p .Tflash) equals the heat production at the asperity (τ c .ν.t c .A asp ) minus the heat buffer due to the interacting water volume (V w .ρ w (c pw .T + L w ) yields: where V asp ¼ A asp ffiffiffiffiffiffiffiffiffiffiffiffiffiffi κ th :π:t c p is the heated solid volume and A asp = π.r 2 is the heated area of the asperity.
Therefore, we computed the flash temperature rise at the contacts at equilibrium following: where water density (ρ w ), specific heat (c pw ) and latent heat (L w ) evolved with pressure and temperature through the thermophysical evolution interpolated from data of NIST [25], shown in Fig. 4c-e. κ th is the rock's thermal diffusivity. Parameter values are given in Supplementary Table. 2.
This idealized model accounted for temperature buffering from asperities by the fluid volume in purely diffusive interaction with it. Supplementary Table 2 presents the values used for calculations that resulted in the flash temperatures of Fig. 4a. The following limitations are noticeable: First, this idealized model does not account for reduction of normal stress due to TP of fluid. Instead, we imposed a constant normal stress with increasing slip in order to observe the theoretical flash temperature reached at the asperities. In our experimental conditions, since the laboratory earthquakes nucleated and arrested spontaneously, we expect slip rates to increase during the dynamic stress drop (from τ 0 to τ min ), accommodating most of the event slip. Then, a deceleration of the slipping zone is expected to occur during the healing phase (from τ min to τ f ) and so a very fast reduction of shear heating is expected until temperatures slightly higher than room temperature. Our model accounted only for the first phase, where the fault slips at constant slip rate. During the rupture arrest phase, the reversibility of the vaporization process should account for fast cooling of the melted asperities and significant reduction of the pressurized volume with an increase in normal stress, arresting fault weakening and increasing fault's strength. Second, further considerations of permeability, porosity and other properties of the rock during mechanical changes induced by seismic slip and the rupture passage are not considered in our calculations 31 . Nevertheless, combining the mechanical results recorded dynamically during earthquake rupture, the observed microstructures and our idealized models bring major insights to the interaction between fault fluids and the weakening mechanisms activated thermally during seismic slip. Note that here vaporization of fault water is reflected in the jump in latent heat, which acts as a heat barrier (Fig.4e), therefore the notion of kinetics of this phase transition is not taken into account in this model. Bulk fault temperature model. We considered a one-dimensional macroscopic fault critically stressed at an initial normal stress (σ n ') with a friction coefficient of 0.7 (ref. 21 ). The fault is sheared at a constant shearing rate v over a thin slip zone of thickness wsz 2,15 where the temperature and fluid pressure increase with shear loading. Note that, in the manuscript, the results presented are for a shear zone thickness of 5 μm. For a finite amount of frictional slip (δ), the generated heat ( q ¼ τ: v wsz 2 ) induces a temperature rise T in the slip zone and then diffuse into the surrounding rock wall. If the bulk fault shear heating phenomenon is drained, we will assume that the fluid pressure in the fault is equal to the initial pore pressure imposed. Conversely, if the conditions are undrained, the generated heat will induce a pore pressure rise (Δpf). Therefore, the thermophysical properties of fault water (ρ w , c pw , their derivatives and η) evolve with pressure and temperature 25 . In this model, no chemical reactions are investigated.
In order to quantify the water volume in the fault, a given fault porosity φ is imposed and so the specific mass capacity of the bulk fault is a function of porosity such that: 15 In this model, we considered the energy and fluid mass conservation equations (in a similar manner as that of ref. 15 ) for the energy and fluid mass conservation in the presence of fluids such that: and where α th ¼ l w ðρ w :c pw Þ À1 is the thermal diffusivity of the fluid and l w is the thermal conductivity of the fluid. λ f and λ n are, respectively, the isobaric thermal expansion coefficients of the fluid volume and of the solid pore space. β f and β n are, respectively, the compressibilities of the fluid volume and the solid pore space. And finally, α hy =k(η.β) −1 is the hydraulic diffusivity of the fault where k is the fault's permeability, η is the fluid's viscosity and β is the compressibility of the fluid volume.
Note that, as discussed by Chen et al. 15 , the latent heat (here L w ) being constant for pressures lower than that of the supercritical phase transition, it vanishes when deriving the energy conservation equation since ∂ρ w :h w ∂T ¼ ∂ρ w ðc w :TþL w Þ ∂T ¼ ρ w :c w , where h w is the enthalpy of water.
In addition, this model does not consider the kinetics of the vaporization transition but rather an instantaneous phase change when temperature is high enough to overcome this transition.
We solve Eqs. (6)-(8) by an explicit finite difference method. On the spatial boundaries, we impose a no-flow condition. The spatial size of the fault and space step (following y axis) is taken such that pressures and temperatures have reached a constant value far from the boundaries at final time.
When fault slip is a purely drained phenomenon, pore fluid pressure is constant in time on the fault during slip. In this case, our model accounts for the temperature rise at a shear stress across the fault taken equal to the shear strength reached in our experiments (τ 0 = 70 MPa). Then the main assumption is that the shear stress follows the effective stress law such that: τ ¼ μ 0 σ n À pf ð Þ , therefore if the deformation is drained, the stress remains constant. This constant shear stress assumption is not fully realistic but it also gives an upper bound to the temperatures reached during fault slip. The results of this model are presented in Fig. 5a and commented in the main text.
When shear heating due to fault slip is a purely undrained-adiabatic phenomenon, pore fluid pressure evolves with time in the fault during slip, and so the thermophysical properties of pore fluid evolve with rising temperature and pressure. In this case, the shear stress along the fault once again follows the effective stress law described before. Therefore, any increase in pore fluid pressure induces a reduction in effective stress. This case gives us a lower bound for the reached temperature rise if TP is the only weakening mechanism activated during slip. The results of this model are presented in Fig. 5b, c and commented in the main text.
Extrapolation to upper crustal depths. In order to extrapolate the model to upper crustal depths, we substituted the experimental stress used in the heat source term of (Eq. 7) for a mean stress taken as the lithostatic overburden gradient of 27 MPa km −1 with a fluid pressure gradient of 10 MPa km −1 , an initial friction of 0.7 and a geothermal gradient of 30°C km −1 (ref. 27 ). The results of this extrapolation are presented in Fig. 6 and commented in the main text.
The following limitations are noticeable: The kinetics of the vaporization reaction are not considered in this model. For details on some attempts to constrain such kinetics, refer to ref. 15 . Instead, we have considered that the vaporization reaction is instantaneous and that the latent heat acts as a heat barrier. In the case of the supercritical transition, the terms concerning the kinetics of the reaction vanish 15,16 . In addition, all heterogeneities that exist in fault zones 1 (in terms of thermal, mechanical and hydraulic properties normal and parallel to the fault plane) are neglected in this model. Data availability. Data are available in supplementary materials. Any further information can be requested from the corresponding author.