Dilatancy toughening of shear cracks and implications for slow rupture propagation

Dilatancy associated with fault slip produces a transient pore pressure drop which increases frictional strength. This effect is analysed in a steadily propagating rupture model that includes frictional weakening, slip-dependent fault dilation and fluid flow. Dilatancy is shown to increase the stress intensity factor required to propagate the rupture tip. With increasing rupture speed, an undrained (strengthened) region develops near the tip and extends beyond the frictionally weakened zone. Away from the undrained region, pore fluid diffusion gradually recharges the fault and strength returns to the drained, weakened value. For sufficiently large rupture dimensions, the dilation-induced strength increase near the tip is equivalent to an increase in toughness that is proportional to the square root of the rupture speed. In general, dilation has the effect of increasing the stress required for rupture growth by decreasing the stress drop along the crack. Thermal pressurisation has the potential to compensate for the dilatant strengthening effect, at the expense of an increased heating rate, which might lead to premature frictional melting. Using reasonable laboratory parameters, the dilatancy-toughening effect leads to rupture dynamics that is quantitatively consistent with the dynamics of observed slow slip events in subduction zones.


Introduction
Dilatancy is a well documented process in intact, healed or overconsolidated rocks: due to shear deformation, microcavities (cracks, grain junctions) open, which induces an overall increase in porosity (Paterson and Wong, 2005, Section 5.3). In fluid saturated rocks, this increase in porosity has the potential to induce fluid pressure drops if the dilatant region is insufficiently drained, producing an increase in effective stress and thus strengthening the material. This phenomenon is called dilatancy hardening.
The importance of dilatancy hardening in fault mechanics has long been recognised in laboratory experiments and in theoretical rupture models. One of the first experimental evidence of dilatancy hardening in crustal rocks was obtained by Brace and Martin (1968) in diabase and granite, who showed that failure strength was increasing with increasing deformation rate as the deformation conditions became more undrained. Similar observations have been made on a range of low-porosity rocks (e.g., Rutter , 1972;Chiu et al., 1983;Duda and Renner , 2013). Martin (1980) further demonstrated that rock failure could be stabilised due to dilatancy, by considerably slowing down the deformation leading to strain localisation and fault slip. New laboratory results by Aben and Brantut (2021) have confirmed the direct stabilisation of ruptures due to shear-induced dilation. Such work was focused on initially intact rock, where dilatancy is occurring in the bulk as well as within the incipient fault zone. In preexisting fault zones, dilatancy can be due to overriding asperities (along bare rock surfaces) and granular rearrangements (if gouge layers are present), leading to net fault zone opening during slip. Such behaviour has been thoroughly documented in artificial fault gouges, for which dilatancy can be related to frictional state evolution (e.g., Marone et al., 1990;Sleep, 2006;Samuelson et al., 2009). Recent laboratory studies have shown that fault zone dilatancy (and compaction) exerts a first-order control on fluid pressure (Lockner and Byerlee, 1994;Faulkner et al., 2018;Brantut, 2020) and slip dynamics (Proctor et al., 2020;Aben and Brantut, 2021).
From a theoretical point of view, the role of dilatancy hardening in the propagation of shear rupture was initially investigated by Rice (1973), who showed that dilatancy tends to increase the fracture energy necessary to drive rupture propagation. Rice and Rudnicki (1979) demonstrated that dilatant hardening leads to a period of stable, quasi-static deformation prior to shear rupture instability. In a one dimensional (spring-slider) fault zone model in the context of slip-dependent strength and dilation, Rudnicki and Chen (1988) showed that dilatancy could prevent unstable fault motion. A similar stabilisation effect was found by Segall and Rice (1995) for faults governed by rate and state dependent friction and dilation. For spatially extended fault slip, the competition between dilation, fluid flow and frictional weakening has been shown to produce slowly propagating slip events, in a manner and parameter range consistent with geophysical observations of slow slip in subduction zones (Liu and Rubin, 2010;Segall et al., 2010;Liu, 2013). Dilatancy has also been shown to promote quasi-static rupture propagation in coupled thermo-hydro-mechanical models where the weakening effect of thermal pressurisation of pore fluid is limited by the strengthening effect of fault zone dilation (e.g., Suzuki and Yamashita, 2007. Most of the aforementioned rupture models including the effect of fault zone dilation have relied on numerical solutions for either the quasi-static (Segall et al., 2010) or dynamic (Suzuki and Yamashita, 2008) equation of motion, effectively solving the fully coupled, nonlinear fracture problem without relying on specific assumptions regarding rupture process zones and weakening behaviour. Liu and Rubin (2010) developed an energy balance approach based on fracture mechanics to determine analytical estimates for slip rate during slow slip transients, which compared well with their numerical results. The success of this approach is further justified by recent work by Barras et al. (2020); Garagash (2021) which demonstrated, in the context of rate and state friction laws, that nonlinear solutions to rupture problems could be recast in the framework of linear elastic fracture mechanics, where rupture propagation is governed by a Griffith-type criterion G = G c , where G is the energy-release rate and G c is a fracture energy that depends on rupture speed and parameters of the nonlinear friction law. Such a connection between nonlinear problems and linear elastic fracture mechanics equivalents was already highlighted by Ampuero and Rubin (2008); Hawthorne and Rubin (2013) for rate and state friction law. Estimates of fracture energy were also given by Segall et al. (2010, §78, 79 and 80) for the coupled, quasi-static, rate and state friction and dilatancy-diffusion problem. While those estimates are very useful to analyse simulation results a posteriori, they are expressed as a function of slip rate, which is a priori unknown in fracture problems. In order to use fracture energy in a Griffith energy balance G = G c , or, equivalently, in a balance between stress intensity factor K and toughness K c , and arrive at elementary predictions for rupture tip dynamics (a time and space history of the crack tip position), fracture energy or toughness would need to be expressed as a function of rupture speed instead of fault slip rate (as done recently by Garagash, 2021).
In addition to its potential impact on stability and fracture energy, one other key aspect of dilatancy in the context of rapid fault slip is its potential to counteract thermal pressurisation of pore fluids and to accelerate heat production. This effect was explored extensively by Garagash and Rudnicki (2003), and some solutions of the coupled thermo-hydro-mechanical problem were exposed by Segall and Bradley (2012). For dilatancy occurring much more rapidly than shear heating, it amounts to resetting the pore pressure inside the fault at the onset of slip (Rice, Figure 1: Schematic of the model setup, with dynamically propagating rupture at speed v r along the x axis under anti-plane geometry. The internal structure of the fault is considered to be a porous material (gouge) of finite width (in the y direction), undergoing shear and dilation.
2006, §50), which increases the potential maximum temperature achieved during slip. This phenomenon plays a major role in determining the onset of frictional melting in low porosity, consolidated rocks (Brantut and Mitchell , 2018); in granite, dilatancy can indeed be so large as to rapidly decrease pore pressure down to vapor pressure during faulting and slip (Brantut, 2020). Recent laboratory measurements of fault zone dilatancy and hydraulic properties during dynamic rupture events (Brantut, 2020) highlight the need for a comprehensive assessment of the effect of dilatancy on rupture dynamics. The present study focuses attention on the effects of fault zone dilatancy on the weakening and thermal response of the fault during dynamic slip, and aims to find a simple description of those effects in terms of fracture energy, to be used in a Griffith-type energy balance. Based on recent laboratory results (Brantut, 2020;Aben and Brantut, 2021), a simple slip-dependent weakening and porosity change model is used (similar to the approach of Rudnicki and Chen (1988)). Although more advanced parameterisations of weakening and dilation based on rateand-state descriptions of friction have been determined experimentally at slow slip rates and used in theoretical investigations (e.g., Marone et al., 1990;Segall and Rice, 1995;Sleep, 2006;Samuelson et al., 2009), a direct dependence of friction and dilation on slip is a limit case of rate-and-state descriptions in response to sudden changes in slip rate, which is particularly relevant to rapid propagation of shear ruptures. A steadily propagating dynamic crack tip model is employed to arrive at semi-analytical results for shear stress, pore pressure and temperature evolution along the propagating rupture. Keeping with the eventual goal of finding a simple, usable form of fracture energy to use in a Griffith energy balance, a particular attention is paid to the contribution of fault zone dilation to fracture energy as a function of rupture propagation speed. This analysis leads to a refined description of the various contributions to fracture energy and energy release rate of different weakening processes, including "near-tip" mechanisms such as intrinsic frictional weakening or dilatancy, which contribute to G c , and "far-tip" mechanisms like pore fluid diffusion or thermal pressurisation (or decomposition, or melting), which are expected to contribute to G.

Elastodynamics
Our main focus is to determine the relative contributions of slip-weakening, dilatancy and thermal pressurisation in the energy balance of dynamic ruptures. A fully dynamic approach would Figure 2: (a) Constitutive laws for friction and porosity evolution with slip (Equations 3 and 4); (b) Propagating crack model, with stress change due to both friction and pore pressure changes (Equation 2). The tip is mostly undrained, and recharge occurs at large distances.
require extensive numerical simulations, imposing somewhat arbitrary choices of initial and boundary conditions. For the purpose of making estimates of fracture energy and dissipation consistent with elastodynamics, it is useful to consider only the tip region of a dynamically propagating, steady-state, semi-infinite shear rupture (Viesca and Garagash, 2015). The steady-state assumption corresponds to an assumption of constant rupture speed v r ; it is valid when rupture speed varies only smoothly as a function of time, over distances (or timescales) much larger than the characteristic dimensions of the non-linear tip region. In a two dimensional configuration where the fault lies along the x axis (Figure 1), the shear stress τ at position X = v r t − x (in a reference frame moving with the rupture tip, see Figure 2b) is related to the slip rate distribution V (X) as whereμ is an apparent elastic modulus that depends on loading mode and rupture velocity. For semi-infinite ruptures, stress drop is neglected compared to strength variations and τ b is set equal to the far-field residual shear strength on the fault. For simplicity, we only consider mode III (anti-plane) loading here; in this caseμ = F (v r /c s )µ where µ is the shear modulus of the fault walls, c s is the shear wave speed, and F (z) = √ 1 − z 2 (e.g., Rice, 1980). In the moving reference frame at constant speed v r , the partial time derivatives are effectively derivatives with respect to coordinate X, ∂/∂t = (1/v r )∂/∂X, so that the dependent variables required to solve the dynamic problem are functions of X only.

Shear strength
In the slipping part of the fault, the shear stress in Equation 1 must equal the shear strength of the fault. We assume that the fault is of finite width, filled with a porous material ( Figure  1). It is initially saturated with a pore fluid, considered chemically inert at the timescale of the phenomena of interest here, at a uniform pressure p = p 0 . The shear resistance to sliding, τ f , is given by Terzaghi's effective stress principle, where σ = σ n − p is the effective normal stress, σ n is the applied normal stress, considered constant, and f is the friction coefficient. The fault zone is assumed to be made of a consolidated material, for instance a healed cataclasite or fault gouge, so that shear strain produces significant (1) "direct" weakening of the fault by a decrease of the friction coefficient from a peak f p to a residual value f r , and (2) dilation of the fault, by way of microcrack opening and overriding of micro-and meso-scale asperities (e.g. Barton, 1976).

Slip weakening behaviour
The former "direct weakening" effect has been thoroughly documented experimentally in consolidated or healed rocks (e.g., Karner et al., 1997;Nakatani and Scholz , 2004), and can be thought of as an irreversible decrease in cohesion; while the physical processes leading to the progressive degradation of strength with ongoing slip remain difficult to constrain in detail, it suffices for the present analysis to capture such a weakening effect as a simple progressive slip-weakening behaviour for the friction coefficient (see a similar analysis in Rudnicki and Chen (1988) where δ denotes the total slip across the fault, and δ w denotes a characteristic slip-weakening distance.
In the context of fast slip (with no subsequent healing), the frictional weakening as phenomenologically described by (3) is not fundamentally different from more elaborate formulations based on rate-and-state friction laws including rapid weakening at high slip rate, for instance due to flash heating (Noda et al., 2009;Brantut and Rice, 2011;Viesca and Garagash, 2015). For dynamic weakening, peak friction is typically of the order of 1 (or, 0.6 to 0.8 if we strictly follow Byerlee's rule of thumb), and weakened friction can be as low as 0.1, for instance due to asperity-scale thermal weakening (e.g., Rice, 2006;Di Toro et al., 2011). The slip weakening distance δ w for dynamic weakening could be as small as a few 10s of µm (e.g., Noda et al., 2009;Brantut and Rice, 2011;Viesca and Garagash, 2015), but there is still considerable uncertainty in choosing this parameter in the phenomenological description (3).
At slower, sub-seismic slip rates, frictional weakening as in (3) is consistent with commonly used rate-and-state friction laws in rate-weakening materials and in the absence of healing: a sudden increase in slip rate does produce an exponential decay of the friction coefficient from a peak (proportional to the log of the velocity increase) to a residual value (proportional to the new "state" of the interface) (Scholz , 2002, , chap. 2).
If we focus our attention to slip-weakening behaviour observed in association with decreases in cohesion (e.g., during failure of intact or consolidated rocks), then appropriate values for peak frictional strength could be as high as 1.5, and the weakened friction coefficient is then of the order of 0.6 to 0.8 (Byerlee, 1978). The slip weakening distance δ w typically associated with the failure process is of the order of 0.1 to 1 mm (e.g., Wong, 1982;Lockner et al., 2001;Ohnaka and Shen, 1999;Ohnaka, 2003;Aben et al., 2019).

Slip-dependent dilatancy
Fault zone dilation due to shear deformation is modelled following the approach of Rudnicki and Chen (1988), and we consider that the inelastic porosity change of the fault material depends directly on slip as ( Figure 2a) where ∆Φ max is the maximum inelastic porosity change at large slip, and δ D is a characteristic slip distance associated with the change in porosity. A phenomenological relationship like (4) is supported by experimental data in consolidated rocks (e.g., Barton, 1976;Teufel , 1981;Brantut, 2020;Aben and Brantut, 2021). In this description, the parameter ∆Φ max can be thought of as the porosity change achieved when the fault zone material reaches critical state, i.e., the stage at which further shear deformation does not impart any bulk volume change. In the case of dilatancy produced by shear rupture in consolidated rocks, the characteristic slip δ D is likely commensurate to the slip weakening distance δ w , since both the weakening and the dilatancy phenomena are driven by the same underlying processes (microcrack opening, linkage and slip). This was the case reported by Brantut (2020); Aben and Brantut (2021) in initially intact granite, where δ D of the order of 1 mm was determined from volume change vs. slip data.
The porosity evolution given by (4) is also compatible with the rate-and-state description given by Segall and Rice (1995) for faults subject to sudden changes in slip rate. In that framework, ∆Φ max is given by a dilatancy factor, of the order of 10 −3 , multiplied to the log of the velocity step, and the characteristic slip δ D is equal to the characteristic slip for state evolution, of the order of 10s of µm.
Overall, the phenomenological description (4) is consistent with both shear failure of consolidated materials and slip across preexisting discontinuities or gouge layers (in the "no-healing" limit), provided that the value of parameters ∆Φ max and δ D are chosen accordingly. Such a versatile description, with a limited number of parameters, is useful to capture the essence of the dilatancy phenomenon.

Fluid flow and shear heating
Due to dilatancy, frictional heating and hydraulic diffusion, the fault zone pore pressure p in (2) evolves as a function of time t during slip. The governing equations for thermal pressurisation have been given numerous times in the literature, and only essential steps are recalled here. All notations directly follow those of Rice (2006); Garagash (2012). The conservation of fluid mass leads to the following governing equation for p (e.g., Rice, 2006): where Λ is the thermal pressurisation coefficient, Θ is temperature, α hy is the hydraulic diffusivity and β * is a hydraulic storage capacity of the fault zone material. The local inelastic change in porosity is denoted n pl . The coordinate y is oriented perpendicular to the fault plane ( Figure  1), and we neglect along-fault fluid diffusion. The temperature evolution of the fault zone is given by (Rice, 2006) whereγ is the shear strain rate, ρc is the heat capacity and α th is the thermal diffusivity. The strain rate is assumed to follow a Gaussian profile of fixed width h across the fault (Garagash, 2012) where V is the slip rate. The governing equation (5) requires a local evolution of porosity at spatial position y, and we assume here that the local porosity rate ∂n pl /∂t is directly proportional to the strain rate, with a proportionality factor given by the derivative of (4) with respect to slip: The formulation (8), together with (7), is consistent with expression (4) for the spatially averaged porosity across the fault. Note that the product ∆Φ max × h gives the maximum total fault zone opening due to inelastic dilation.
The solution for pore pressure at y = 0 as a function of time is given in integral form as (Garagash, 2012) where ∆Φ is the derivative (4) with respect to slip. The kernels K and A are given by Garagash (2012, , Appendix A). The coupling between pore pressure and temperature introduces two diffusion timescales, a thermo-hydraulic one, , and a purely hydraulic one, T hy = h 2 /4α hy . In addition, thermal pressurisation introduces a natural weakening length scale which can be chosen as δ c = ρch/(f r Λ). Here, some flexibility exists in the definition of δ c since friction coefficient is not constant; the choice of f r is made to model the case when thermal pressurisation becomes significant only at slip distances greater than δ w , beyond which friction is in its weakened state (a similar assumption was used by Noda et al., 2009;Viesca and Garagash, 2015, in the case of early weakening by flash heating).

Dilatancy toughening of slip-weakening faults
Let us first consider the rupture problem without thermal pressurisation, by assuming that Λ is very small, so that the slip weakening distance δ c is much larger than both δ w and δ D . Such a situation would arise in fault zones with large pore space compressibility, which seems to be the case of fresh granitic faults as reported by Brantut (2020).

Effect of dilatancy on crack tip stress and slip rate
The coupled problem (2), (9) and (1) can be reformulated in terms of normalised quantities x/L w , ∆τ /τ r , V /V w and δ/δ w , where ∆τ = τ − τ r is the stress drop and The initial effective stress is denoted σ 0 = σ n − p 0 . There are 4 parameters controlling the behaviour of the system, the ratios where we denote T = L w /v r . The ratio T hy /T = v r T hy /L w compares the hydraulic diffusion time to the characteristic time over which the slip-weakening region propagates along the fault, and is the key parameter controlling the transition from drained (small T hy /T ) to undrained (large T hy /T ) conditions near the rupture tip. The solution is determined using the quadrature method detailed in Viesca and Garagash (2018). A set of solutions is presented in Figure 3, where constant δ D /δ w = 1 and f p /f r = 2 were chosen as representative values. For significant values of ∆Φ/β * σ 0 (i.e., a nonnegligible undrained de-pressurisation due to dilatancy), we observe a clear transition from a drained regime at T hy /T 1, where the crack tip behaviour is essentially governed by the dry slip weakening mechanisms, to a dual scale behaviour at T hy /T 1 where the tip is undrained, governed by a slip-dependent behaviour (combination of slip weakening and undrained effective stress law), and an extended region (at x > v r T hy ) governed by the progressive pore pressure recharge of the fault zone from the off fault regions (see also Figure 2b). This second weakening 0.0 0.5 where k is constant that depends on the undrained strength and is obtained numerically for each solution; Blue dashed line is the LEFM asymptote V (x)/V w = 2/ √ π × (x/L w ) −1/2 for the drained (dry) tip behaviour; Orange dashed lines correspond to the purely diffusive asymptote given in Equations (12) and (13). behaviour is purely diffusive, and at x v r T hy the strength is well approximated by (this follows from approximation of the second integral in (9), assuming concentrated source at t = 0): and the slip rate behaves as (this follows from inversion of the one-sided Hilbert transform in (1) and direct computation of the resulting integral) With increasing undrained pore pressure change (parameter ∆Φ max /(β * σ 0 )), the peak slip rate at the crack tip decreases, due to a smaller net strength drop under undrained conditions, while the slip rate in the extended "far-tip" region controlled by diffusive pore pressure recharge increases.

Stress intensity factor and fracture energy: Small scale yielding approximation
Now equipped with a numerical solution and relevant asymptotes for the stress distribution along the moving crack, we can turn our attention to the desired analysis in terms of integrated fracture mechanics quantities. By construction, there is no stress or slip rate singularity at the tip of our moving crack. However, we can consider the situation where all the strength change due to the coupled slip weakening and dilation process occurs over a length scale much smaller than the finite size of the expanding crack. In that case, the stress distribution corresponds to a stress intensity factor that needs to be overcome for the crack to keep expanding, which is given by (Kostrov , 1966) In the quasi-static case (v r c s ), it simplifies to (Rice, 1980, pp. 598, 604) The purely drained end-member case, where ∆τ is an exponentially decreasing function of slip, provides a lower bound K drained c for the stress intensity factor. This lower bound is obtained from the finite fracture energy associated with the exponential slip-weakening process (G drained where the relationship between stress intensity factor K and energy release rate G was used. When dilatancy produces a significant pore pressure drop, the stress in the vicinity of the rupture tip increases, and thus K c increases. The diffusive nature of the stress evolution introduces some surprising complexity in the problem: the integral of (15) does not converge. Following Rice (1973), a cutoff length cut is introduced beyond which the strength is considered to be exactly zero, in accordance with the original assumption that all the strength evolution remains confined to a small region near the crack tip. The stress intensity factor computed this way exhibits a marked increase with increasing diffusion time relative to frictional weakening time T hy /T , that is, with increasing rupture speed ( Figure 4). Therefore, dilatancy produces an effective toughening of the fault, which potentially limits the rupture speeds. Using the approximation (12) for the stress distribution near the crack tip, the dynamic stress intensity factor is approximated by This expression for stress intensity factor has a similar behaviour to that derived by Rice (1973) for the case of an infinitesimally thin slipping zone and quasi-static rupture growth, where K was found to increase as √ v r , and a weak dependence on cut was established.
A complementary approach to computing a critical stress intensity factor is to estimate the equivalent fracture energy associated with the weakening processes at the crack tip. Following Rice (2006), we can use a generalised form of fracture energy, better called "breakdown work", defined as which reduces to the conventional form established by Palmer and Rice (1973) when a well defined constant residual stress exists and slip is large compared to any characteristic slipweakening distance. Under drained conditions, for instance at very small rupture speed (T hy /T 10 3 10 2 10 1 10 0 10 1 10 2 10 3 Figure 4: Stress intensity factor required to drive crack expansion, normalised by its perfectly drained value, as a function of hydraulic diffusion time relative to characteristic frictional weakening time. The cutoff length was set to cut = max{100L w , 10 × v r T hy } in the numerical computation. T hy /T increases with increasing rupture speed. Blue dashed line indicates the purely drained limit. Orange dashed lines show the analytical estimate as computed from the sum of the drained end-member and the large T hy /T end-member (18), using cut = 10 × v r T hy .
1), G becomes almost constant at slip distances that are large compared to δ w ( Figure 5). In that case, G ≈ G drained c and there is a well defined cohesive zone of small dimensions. With decreasing drainage (increasing T hy /T ), G evolves in two distinct regimes: (1) an increase and stabilisation up to slip distances of the order of δ w , and (2) a marked increase with increasing slip beyond δ w . The first regime corresponds to a purely undrained behaviour, where shear strength evolves as The undrained behaviour at the tip is associated to a well defined, constant fracture energy and the subsequent increase occurs due to diffusive recharge of the fault zone. The breakdown work, as computed by (19), does not converge to a fixed value at large slip: this is consistent with the theoretically infinite toughness and the requirement of a cutoff length cut that limits the size of the cohesive zone. The unlimited increase in G exhibited by our slip-weakening, dilatant fault zone model, is a general feature of weakening models limited by hydro-thermal diffusion processes (Rice, 2006;Viesca and Garagash, 2015;Brantut and Viesca, 2017). For rupture propagating in consolidated rocks, using representative parameter values of δ w ≈ δ D ≈ 1 mm, µ ≈ 30 GPa, (f p /f r − 1)τ r ≈ 30 MPa and α ≈ 10 −5 m 2 s −1 , the key controlling parameter T hy /T ranges from around 10 −4 at v r ∼ 1 km/day up to 10 at v r ∼ 2 km/s when considering a relatively thin dilatant zone of h = 1 mm, and from 1 to 10 5 over the same rupture speed range when considering a thicker zone of h = 10 cm. For v r approaching c s ,μ vanishes and thus T hy /T becomes unbounded. Therefore, for thin shear zones, dilatancy is expected to produce dramatic toughening only at dynamic rupture speeds, whereas large toughening is predicted even at small rupture speeds, of the order of 10 km/day, if dilatancy occurs in thick zones. 3.3 Stress intensity factor and fracture energy: "near" and "far"-tip contributions While the toughening effect is somewhat intuitive, one caveat is that the concept of toughness and the relevance of the stress intensity factor to characterise crack expansion may be limited, since it relies on the small scale yielding assumption: in computing K c and using it in terms of "crack resistance", one has to assume that the nonlinear end region, where the stress evolves from the nominal strength to the residual value, is small compared to the overall crack size. The need for a cutoff length cut to prevent the mathematical divergence of the integral in (15) indicates that the nonlinear end region may not be small in general: it is only approximately true when the rupture dimensions are much larger than the drainage length scale (v r T hy ). In practice, for finite ruptures, the integral would extend along the full length of the rupture, making the toughness size-dependent (with longer ruptures having larger toughness than short ones). Size-dependent toughness is consistent with the results depicted in Figure 5 showing slip-dependent breakdown work: longer ruptures, which typically correspond to larger slip, have larger breakdown work. Such a size-dependence of toughness is commonly observed in engineering materials due to so-called crack bridging effects (e.g., Cox and Marshall , 1994). Here the matter is further complicated by the time-dependence of the strength evolution, which would make the toughness time-dependent as well.
An alternative to grouping all nonlinearities in the toughness and making it size-and timedependent, is to separate the various contributions to stress intensity factor occurring at different scales. If some of the strength drop is not confined to a small region near the crack tip, which is the case when the rupture dimension is comparable to or smaller than the diffusion length scale, it will not contribute to toughness (or fracture energy) per se, but to the stress intensity factor (or energy release rate). Here, it is natural to have the nonlocal contribution of diffusive fault zone recharge into G, and keep the undrained contribution of dilatancy into G c . In any case, the energy balance G = G c is always verified, but G now includes nonlinearities due to the constitutive behaviour of the fault.
A simple illustrative example can be worked out: consider a semi-infinite crack along direction x, propagating at constant rupture speed v r in a background stress field given by τ b > 0 for x ≥ 0 and τ b = 0 for X < 0. The stress drop along the crack is τ b − τ f , and for sufficiently large diffusion time, the strength evolution τ f can be split into a near tip, undrained contribution, and a spatially extended diffusive contribution. The near tip undrained contribution corresponds to an undrained toughness K undrained c = 2μG undrained c . This toughness must be matched by the stress intensity factor associated with the background loading, offset by the spatially extended diffusive contribution (12). This stress intensity factor, computed from (14), yields In the limit v r T hy , the stress intensity factor is simply: so that dilatancy has the effect of reducing the stress drop, and limiting the available energy for crack tip extension. By contrast, at large rupture dimensions ( v r T hy ), the stress intensity factor becomes and the contribution of dilatancy eventually becomes negligible compared to that of the applied load.
The simple example of a semi-infinite crack propagating at constant speed is clearly not a fully realistic earthquake or slow-slip model, and the choice of applied load was made for mere mathematical convenience. However, it is quite illuminating to establish that dilatancy has a two-fold effect: (1) changing the fracture energy at the crack tip, and (2) reducing the stress drop and therefore the energy release rate in a time-dependent, diffusive manner. The different contributions of dilatancy and drainage can be recast in the Griffith energy balance, which provides an equation of motion for the crack tip through the dependency of G on v r and crack tip position. The energy release rate G is thus impacted by the combination of applied load and diffusive pore pressure recharge along the crack, while the undrained behaviour (combination of slip weakening and pore pressure drop) controls fracture energy. The energy release rate depends quadratically on K, and is therefore not simply offset by the diffusive pore pressure recharge far from the tip. Thus, in general, the full contribution of dilatancy and drainage to the dynamics of rupture is not simply a change in fracture energy. The approximation of small scale yielding made in the previous section (see Equation 18), which followed the developments of Rice (1973), is thus valid only in the limit of ruptures that extend well beyond the region impacted by pore pressure recharge.

Comparison to dynamic simulations
Whether due to a direct toughening effect, or to a reduction in energy release rate, dilatancy is expected to slow down rupture propagation. This can be tested in full elastodynamic simulations (see details in Appendix A). An illustrative example is shown in Figure 6, where key constitutive parameters were set to ∆Φ max /(σ 0 β * ) = 0.5, f p /f r = 2, and the background stress outside the nucleation zone was set to 0.5τ r above the residual frictional strength. In that case, a purely undrained behaviour would amount to zero stress drop, so that rupture should not grow unless drainage becomes significant (G, or K, is exactly zero under purely undrained conditions, cf.  Figure 6: Slip rate (a) and stress (b) profiles during an elastodynamic rupture simulation with ∆Φ max /(σ 0 β * ) = 0.5, δ w = δ D , f p = f r , and a normalised hydraulic diffusion time T hy /T = 10 3 , where the characteristic frictional weakening time is defined as T = L w /c s , the characteristic distance is L w = δ w τ r /µ, and the characteristic slip rate is V w = δ w /T . The background stress was uniform at 1.5τ r , except for a small region of size 12L w near the origin where its value was 1% above the peak strength. The rupture is symmetric with respect to x = 0 and only the domain x > 0 is displayed. Profiles are plotted at time intervals equal to 12.5T . Blue lines in each plot show the prediction of rupture tip trajectory from the Griffith energy balance.

Equation 23
). In the simulation, the rupture tip is indeed arrested at early times, at around x/L w = 18 (where L w = µδ w /τ r ). The only mechanism that allows rupture to grow any further is the diffusive pore pressure recharge of the fault zone, which increases the strength drop in the interior of the crack, and increases the energy release rate. After its momentary arrest, rupture slowly accelerates, but remains well below the shear wave speed (the limiting speed in mode III) over the full simulation time, of the order of 400 × µδ w /(τ r c s ). By comparison, simulated ruptures in the same background stress but without dilatancy accelerate to v r ≈ c s almost immediately after nucleation.
The stress and slip rate near the tip of the spontaneous dynamic rupture are well characterised by the steady-state solution using the appropriate transient rupture speed ( Figure  7). The tip is indeed under essentially undrained conditions, and the slip rate is reasonably well characterised by the linear elastic fracture mechanics (LEFM) solution V /V w = 2K undrained c / √ π × (x/L w ) −1/2 , corresponding to the undrained stress intensity factor as defined by K undrained c = 2μG undrained c . In the dynamic model, the progressive strength drop associated with drainage does not lead to a significant increase in slip rate further away from the tip (as predicted by the semi-infinite crack solution) due to the finite crack size of the model.
The usefulness of the LEFM approximation can be tested by comparing the rupture tip trajectory obtained from the full numerical solution to that based on the crack tip equation of motion arising from the Griffith energy balance. In view of the undrained behaviour at the tip, it is appropriate to equate the dynamic energy release rate to the undrained fracture energy G undrained c , and to include the strength change due to pore fluid diffusion (given asymptotically by Equation (12)

Link with thermal pressurisation
In the previous Section we have established how dilatancy and pore fluid diffusion influence the dynamics of rupture, without regards to thermal effects. In general, thermal pressurisation is expected to produce significant weakening for slip distances larger than δ c . In granite gouge, δ c is typically of the order of 10 mm (Brantut and Platt, 2017), which is around one order of magnitude larger than δ w or δ D . In freshly fractured granite, Brantut (2020) reported very high pore space compressibility compared to previously reported "mature" fault gouge material, which indicates that δ c is likely much larger than the estimate of 10 mm in consolidated rocks.
In any case, thermal pressurisation is not just superimposed onto the slip weakening and dilatancy processes, but is coupled with them: as noticed by Garagash and Rudnicki (2003), the de-pressurisation effect of dilatancy could be compensated by increased shear heating and fluid pressurisation.
When dealing with the fully coupled slip weakening, dilatancy and thermal pressurisation dynamic problem, a total of 6 nondimensional quantities control the system's behaviour: Let us consider a representative case by choosing δ D /δ w = 1, f p /f r = 1, χ = 10, δ c /δ w = 10, and examine how dilatancy and hydro-thermal diffusion influence rupture dynamics, solving the steady-state problem (Equation (1) with (9), (3) and (4)). When thermal pressurisation is active, the residual strength is zero, so that the correct assumption to solve the semi-infinite rupture problem is that the applied stress is also zero.
In the absence of dilatancy and with relatively fast diffusion (T hy /T = 1; Figure 8(a)), the stress evolution follows two stages of weakening, one associated with frictional weakening down to around τ /τ r = 1 over a distance of the order of L w , and a subsequent, slow decrease down to near 0 strength, associated with thermal pressurisation. This is the situation presented in detail by Viesca and Garagash (2015, their figure 1b). When dilatancy is included, here choosing ∆Φ max /(β * σ 0 ) = 0.5, a significant strengthening is observed and the first stage of weakening is effectively delayed. However, at large distance from the tip, strength, slip rate and temperature become indistinguishable from the nondilatant case. The weakening due to thermal pressurisation is indeed similar between the two cases due to (1) the extra heating rate (temperature is initially higher in the dilatant case) in the undrained regime, and (2) the diffusive pore pressure recharge from the off-fault region.
The same general pattern is observed when diffusion time is larger (T hy /T = 10 3 , Figure  8(b)): an initial delayed weakening due to frictional slip weakening coupled to dilatancy, followed by weakening driven by thermal pressurisation with an increased heating rate compared to the nondilatant case. The effect of dilatancy vanishes at distances x v r T hy . When δ D δ c (as in the case shown in Figure 8), the adiabatic, undrained behaviour including dilatancy is well described by a simple resetting of the initial pore pressure by the quantity ∆Φ max /β * (Rice, 2006, §43); this is illustrated by the temperature rise, which is offset by ∆Φ max /(β * Λ) compared to the nondilatant case (see Appendix C.1), in the domain L w δ c /δ w < x < v r T hy .
The main effect of dilatancy when coupled to thermal pressurisation is a relative strengthening near the crack tip. This near-tip strengthening corresponds to an increase in breakdown work at small slip distances (Figure 9). Under purely undrained, adiabatic conditions (near the crack tip), there is a well defined fracture energy (the integral (19) converges). For constant friction coefficient, the fracture energy is simply where it is easily seen how dilatancy is merely offsetting the effective stress. The fracture energy (27) is a good approximation for the breakdown work at small slip when diffusivity is low (Figure 9, dashed lines). For large rupture dimensions and large slip distances, hydrothermal diffusion becomes significant and the effect of dilatancy on strength and breakdown work becomes progressively negligible (see curves merging in Figures 8 and 9).
Similarly to the situation where thermal pressurisation was neglected, the continuous increase in breakdown work with increasing slip, at large distances from the crack tip, indicates that the weakening may not be confined to a well-defined small region near the crack tip. Thus, the breakdown work should not be given the meaning of "fracture energy" in the sense that it  Figure 9: Breakdown work, normalised by its finite drained value, as a function of slip, for dynamic semi-infinite crack models computed using f p /f r = 2, δ D /δ w = 1, δ c /δ w = 10, χ = 10, T hy /T = 1 or 10 3 , and with (black lines, ∆Φ max /(β * σ 0 ) = 0.5) or without (green lines, ∆Φ max /(β * σ 0 ) = 0) dilatancy. The dashed lines correspond to the fracture energy computed using (27) with and without dilatancy.
may not be used directly in a crack tip equation of motion of the kind G = G c . Rather, as discussed in the previous section, if the weakening scales associated with friction, dilatancy and undrained, adiabatic thermal pressurisation are much smaller than the length scale associated with hydro-thermal diffusion, one may separate the tip behaviour and use the small scale yielding approximation to determine a constant fracture energy (e.g., G adiab.,undr. c ), and the continued weakening far from the crack tip will then contribute to the stress intensity factor (or energy release rate).
Because dilatancy is expected to occur primarily over the first few millimetres of slip, concomitant with any frictional weakening process, its contribution is significant in the near-tip region. The key effect of dilatancy is therefore to increase the undrained fracture energy, which has the potential to locally slow down or stop rupture propagation.

Onset of vaporisation and melting
Throughout this paper the properties of the pore fluid have been assumed constant. This is a reasonable hypothesis when the fluid is a liquid, but the strong changes in pressure (and temperature) associated with dilatancy and thermal pressurisation can significantly impact the fluid compressibility, thermal expansivity and heat capacity. Such changes have been invoked to explain slip dynamics in water saturated granite (Acosta et al., 2018), and recent experimental data have brought direct evidence for fluid vaporisation (with dramatic increase in fluid compressibility) due to dilatancy during rock failure (Brantut, 2020;Aben and Brantut, 2021).
Is pore fluid vaporisation a widespread phenomenon during rupture? A simple assessment was made by Brantut (2020) based on laboratory measurements of stick slip in granite, assuming isothermal and undrained conditions. Vaporisation was predicted to occur at modest slip distances (less than 1 cm) under upper crustal conditions. This simple assessment can be revised by considering fluid diffusion and thermal pressurisation, which tend to limit the pressure drop associated to dilatancy. The maximum pore pressure drop is given by ∆Φ max /β * , but a decreasing fraction of this maximum is achieved with increasing drainage and increasing contribution of thermal pressurisation (Figure 10). In the case of freshly fracture granite reported by Brantut (2020), the slip weakening distance associated to thermal pressurisation is most likely very large: the thermal pressurisation factor Λ is small due to the large pore space compressibility of the newly formed gouge material, of the order of 10 −8 Pa −1 . Combining this large compressibility with typical parameters for granite, using a fault thickness of 3 mm, a friction coefficient of 0.6, a heat capacity of 2.6 MPa/K and a thermal pressurisation factor of Λ = 0.03 MPa/K, we find δ c ≈ 0.4 m. This is indeed much larger than the inferred frictional slip weakening distance and characteristic slip associated with dilatancy, which are of the order of 1 mm. In that case, neglecting thermal pressurisation to compute pore pressure drop is justified, and the maximum pore pressure drop should remain close to the undrained limit ∆Φ max /β * if rupture is sufficiently fast (undrained tip). In laboratory experiments, direct measurements of pore pressure drop during faulting of intact granite are of the order of 30 MPa for slip of the order of 1 mm, so that vaporisation is expected in roughly the upper 3 km of the crust in such rock type (Brantut, 2020;Aben and Brantut, 2021). With increasing temperature, the vaporisation pressure increases up to the critical point, around 20 MPa at 370 • C, so that vaporisation could be expected down to 5 km depth. Such estimates are necessarily coarse, and strong local variations are expected due the natural roughness of faults (e.g., in large pull-apart regions).
When thermal pressurisation is not efficient, as in the case of newly formed faults in granite, frictional heating is large, which can lead to melting of the fault material. The characteristic temperature rise for undrained, adiabatic behaviour in the presence of rapid dilatancy at the onset of slip is given by (σ 0 + ∆Φ max /β * )/Λ (see Appendix C.1). At very large slip distances compared to δ c , the effect of dilatancy becomes negligible (see Appendix C.2). A thorough examination of the onset of melting during thermal pressurisation was given by Rempel and Rice (2006); their analysis remains valid at large slip, and only needs to be amended at small slip by accounting for an initial offset in pore pressure equal to ∆Φ max /β * .
In the case of granite, the slip at which melting occurs is considerably shortened by dilatancy. For instance, at a depth of 5 km, assuming lithostatic normal stress and hydrostatic initial pore pressure, the initial effective stress is σ 0 = 90 MPa. With an average thermal pressurisation factor of Λ = 0.1 MPa/K, intermediate between the low end-member value of 0.03 reported in Brantut (2020) for fresh fractures and a high end-member value of 0.3 MPa/K reported in Brantut and Platt (2017) for mature granite gouge, the adiabatic, undrained temperature rise is of 900 K. Melting is therefore expected only at slip distances several times larger than δ c , (here a few centimetres if we assume a shear zone width of a few millimetres), well into the regime where strength is already significantly reduced by thermal pressurisation. With a pressure drop due to dilatancy of the order of 30 MPa, commensurate to that measured in the laboratory, the adiabatic, undrained temperature rise becomes 1200 K, so that only a fraction of δ c is required to trigger bulk melting of the shear zone, i.e., before any significant strength reduction due to thermal pressurisation.
Thus, dilatancy dramatically facilitates frictional melting. The common occurrence of pseudotachylytes in low-porosity, consolidated rocks (Sibson and Toy, 2006) is therefore expected not only because of the potential "dryness" of the initial material, but because dilatancy would either vaporise any preexisting fluid or trigger fast melting due to the induced pore pressure drop (see also Brantut and Mitchell , 2018, where this is discussed with an example from the field). Parameter values are f p /f r = 2, δ D /δ w = 1, χ = 10, ∆Φ max /(β * σ 0 ) = 0.5, and thermal pressurisation slip weakening distance is changed as reported on the graph.

Diffusion-driven rupture nucleation and growth
As observed in dynamic simulations (see illustrative example in Figure 6), dilatancy has a dramatic impact on the dynamics of rupture during the nucleation phase. A thorough analysis of the transition from aseismic to seismic slip in slip-weakening, dilatant faults has been conducted by Ciardo and Lecampion (2019) in the context of injection-induced slip. Their analysis only considered along-fault fluid diffusion, assuming impermeable fault walls. Here, we considered the alternative hypothesis whereby only across-fault fluid flow was accounted for, which is reasonable when the pore pressure gradient along the x direction (which scales with 1/L w ) is much smaller than that along the fault-perpendicular direction (which scales with 1/h). Despite this difference, one important result of Ciardo and Lecampion's analysis remains valid, and follows from the stress intensity factor computation given in Equation (22): there exists a minimum dilatancy amplitude above which rupture stabilisation can occur, given by At ∆Φ max > ∆Φ c max , the undrained (isothermal) residual strength is above the background stress and rupture cannot propagate unless some fluid flows back into the fault zone. The rate of fault growth is then entirely determined by the rate of pore pressure recharge. The dynamic simulation shown in Figure 6 corresponds to the case ∆Φ max = ∆Φ c max , and a exhibits a large nucleation timescale compared to typical elastodynamic timescales. In addition to the stabilising effect due to offsets in residual strength in undrained conditions, dilatancy also has the potential to produce transient slip-strengthening behaviour at the crack tip, which considerably increases the nucleation size of ruptures loaded with relatively uniform stress (Brantut and Viesca, 2015).
While a thorough analysis considering two-dimensional hydraulic diffusion and finite crack dimensions should eventually be conducted to make quantitative estimates of critical nucleation dimension and time, it is clear that dilatancy has a first order impact on fault stability and rupture initiation. Low porosity rocks, either intact, healed or sealed (e.g., granites) are the most affected by this process.

Dynamics of slow ruptures driven by dilatancy
In models based on rate and state friction laws, dilatancy hardening has been shown to be a viable mechanism for the generation and propagation of slow slip events in subduction zones (Segall et al., 2010). There are many characteristics of slow slip events that need to be explained quantitatively by any proposed physical model: spatio-temporal dynamics of events, recurrence rate, moment-duration relationship, relationship with earthquake location and timing, etc (e.g., Bürgman, 2018). Here, our very simple rupture model is not aimed at explaining all those observations; however, it is useful to test whether dilatancy hardening alone can explain the slowness of rupture propagation during slow slip events, and what parameter range would be required to do so.
Slow slip events in the Cascadia subduction zone typically propagate at speeds of the order of 10 km/day. Their spatial extent ranges from 50 to 100 km, and they are associated to stress drops ranging from 2 to 100 kPa (Schmidt and Gao, 2010;Michel et al., 2019). The stress intensity factor at the tip of such ruptures can be approximated by Hawthorne and Rubin (2013); Weng and Ampuero (2019) where ψ is a geometrical constant (≈ 0.87 in the geometry used by Hawthorne and Rubin (2013), and given by ≈ 2/π for deeply buried antiplane faults according to (Weng and Ampuero, 2019)), ∆τ is the stress drop and W is the width of the ruptured region. Using the spatial dimension and stress drop estimates from geodetic measurements, we obtain K ranging from 0.5 to 24 MPa m 1/2 . This stress intensity factor has to be matched by the toughness that arises from frictional strength. Under purely drained conditions, the toughness is (Equation For slow slip events in subduction zones, we can reasonably assume that the fault zone material is not overconsolidated or well cemented, so that conventional rate and state friction is an appropriate description of the strength. Our simplified model of purely slip-weakening friction corresponds to the "no-healing" end-member of the slip law (Garagash, 2021); in that approximation, the friction drop at the rupture tip is therefore of the order of f p − f r ∼ 10b (Hawthorne and Rubin, 2013;Garagash, 2021), where b is the evolution parameter of the friction law, itself of the order of 0.01. The slip weakening distance is of the order of 1 to 100 µm. The effective stress in the regions experiencing slow slip events is low, of the order of a few megapascals, as inferred by seismological observations (Bürgman, 2018, Section 2.3). Using a shear modulus of µ = 30 GPa, we arrive at K drained c ranging from 0.1 to 1.7 MPa m 1/2 , which overlaps with the K estimate only when very low stress drops are considered in conjunction with large slip weakening distances and effective stress. In addition, rupture events are not expected to be stable if we consider only a balance between K and K drained c , since any increase in K would lead to dynamic acceleration of the rupture tip (Weng, 2021). When dilatancy is accounted for, ruptures are stabilised since any increase in K and rupture speed can be matched by a corresponding increase in K c (Figure 4; see below).
As rule of thumb, for any stabilisation effect to be significant, i.e., for K c to increase significantly above K drained c , the magnitude of the undrained pore pressure drop has to be (at least) comparable to that of the stress drop. If dilatancy is the stabilising process of slow slip events, we therefore expect the associated pore pressure drop to be, at the minimum, in the range 1 to 100 kPa.
More precisely, the contribution of dilatancy and subsequent pore fluid diffusion can be estimated analytically using Equation (18), assuming that the hydraulic diffusion time T hy is large compared to the elastodynamic timescale T = µδ w /v r (f p − f r )σ 0 , that is, when the length scale v r T hy is much larger than the frictional cohesive zone size. Choosing cut = 10 × v r T hy (or, alternatively, cut = W , with little quantitative impact on the result due to the weak  Figure 11: Combinations of dilatancy and diffusion parameters required to satisfy K = K c for slow slip events for a range of stress drops. The frictional contribution to toughness is neglected compared to that of dilatancy. Dashed lines correspond to hydraulic diffusion times such that the characterstic distance v r T hy is smaller than typical frictional cohesive zone size, rendering the use of Equation (18) imprecise.
dependency of K c on the cuttoff length), we can look for the combination of parameters T hy (hydraulic diffusion time) and f r ∆Φ max /β * (dilatant strength increase) required to match the stress intensity factor K estimated from geophysical data. For stress drops of the order of 30 kPa, rupture speeds of 10 km/day can be explained by a combination of f r ∆Φ max /β * = 0.3 to 0.1 MPa and T hy = h 2 /4α hy = 300 to 3000 s, respectively ( Figure 11). This corresponds to an undrained pore pressure drop at the rupture tip of the order of 0.1 MPa, so that only very modest dilatancy is required to severely limit the rupture speed, provided that diffusivity is low enough. For a compressibility β * of the order of 10 −10 to 10 −9 Pa −1 , a pressure drop of 0.1 MPa corresponds to a porosity change ∆Φ max from 10 −3 to 10 −2 %, which is commensurate with the dilatancy parameter inferred by Segall and Rice (1995) from Marone et al. (1990) data, andused by Segall et al. (2010) in their model for slow slip events. Even though the model used here is only slip-dependent and not directly rate-and-state dependent, both approaches point to relatively small dilatant effects (i.e., not of the order of tens of MPa as could be expected when fracturing intact material). Regarding pore fluid diffusion, in clay gouge, Faulkner et al. (2018) report diffusivities of the order of 10 −8 m 2 /s, so that shear zones of around 6 mm width would be sufficiently thick to obtain diffusion time consistent with the model. If diffusivity of the fault zone material is comparable to that reported for exhumed continental faults, α ≈ 10 −6 m 2 /s (Wibberley and Shimamoto, 2003), considerably thicker shear zones of 6 cm would be needed to explain the rupture dynamics.
The parameter estimates obtained from our simple stress intensity factor analysis of slow slip events are realistic and consistent with laboratory measurements. As explained above, the rupture model used here is not designed to reproduce all the characteristics of slow slip events, but it does explain how slow rupture speeds can be achieved with dilatancy and pore fluid diffusion as the main control of strength.
Compared to the rate-and-state dependent approach used by Segall et al. (2010), the slipdependent friction and porosity model we have employed is considerably simpler, but leads to similar rupture dynamics. The computation of K relies on the ruptures being confined in a strip of finite width (Hawthorne and Rubin, 2013), which is not part of the assumptions used by Segall et al. (2010). In a strictly two-dimensional configuration, the stress intensity factor is proportional to the squareroot of the rupture length. To see how ruptures might grow in this geometry, let us consider again, as in Section 3.3, a semi-infinite fault loaded by a background stress τ b > 0 in the region x > 0, and τ b = 0 elsewhere. At a given time t, the rupture has grown to a distance L in the loaded region, so that the stress intensity factor at the tip is In the small scale yielding approximation with cutoff length c = 10 × v r T hy , this stress intensity factor must equate the toughness given by Equation (18), which is proportional to √ v r = dL/dt. Neglecting dynamic effects at small v r , we arrive at the following differential equation for L(t): Assuming, as before, c = 10×v r T hy (keeping in mind that the dependency on cut is very weak), the last term in (31) is a constant, approximately equal to 0.17. The solution of (31) is thus an exponential growth: where the growth rate is Using the parameter values discussed above, with a stress drop of the order of τ b −τ r = 0.01 MPa, a dilatant strength contribution of the order of f r ∆Φ max /β * = 0.1 MPa, and a diffusion time of T hy = 10 3 s, we obtain a growth rate of t g ∼ 5.9 × 10 5 s, i.e., about one week. Such a timescale is commensurate to the total duration of slow slip events (e.g., Michel et al., 2019), and is much longer than any timescale associated with wave propagation along the fault. While the exponential growth might be viewed as unstable, the small growth rate, essentially governed by the diffusion timescale, makes these ruptures far less sensitive to variations in stress (or, energy release rate) compared to ruptures driven by constant toughness (or, fracture energy). In that sense, dilatancy produces relatively stable ruptures, which is consistent with observations of slow slip events in subduction zones.

Conclusions
By analysing a crack model incorporating slip weakening friction and slip-dependent dilatancy coupled to fluid flow, we have established that dilatancy tends to limit crack propagation by increasing the stress intensity factor required for crack growth. At sufficiently large rupture speeds compared to the rate of fluid flow, the rupture tip is undrained, and the frictional strength is larger than in the nondilatant case. Further away from the tip fluid flow (i.e., pore pressure recharge) leads to a progressive decrease in strength back to the drained residual strength.
In a small-scale yielding approximation, the effect of dilatancy amounts to an increase in toughness that scales with the squareroot of rupture speed. This approximation is valid when the rupture dimension is much larger than the along-fault hydraulic diffusion length scale. For shorter rupture dimensions, the strength variations due to pore fluid diffusion are not confined in a small near-tip region, and the approximation of small-scale yielding becomes less accurate. Nevertheless, dilatancy still has the effect of increasing the applied stress intensity factor necessary for crack growth, which can be decomposed in a near-tip contribution (a relative decrease in toughness compared to the nondilatant regime) and a nonlocal contribution (a reduction in the net stress intensity factor) which includes a transition from undrained strength near the tip to the drained strength away from the tip. This natural decomposition of stress intensity factor into identifiable contributions does not have a simple parallel in terms of energy release rate and fracture energy, since G depends quadratically on K, which couples near-tip and "far"-tip effects.
When effective, thermal pressurisation is seen to compensate somewhat the strengthening due to dilatancy through an increase in heating rate. The resulting temperature increase is transiently larger than in nondilatant cases, which can facilitate frictional melting. However, for large slip distances, the effect of dilation on temperature and strength vanishes due to thermohydraulic diffusion with the surrounding medium.
The dilatancy toughening effects analysed here have direct applications in the study of the dynamics of slow slip events. A stress intensity factor analysis shows that the slowness of ruptures propagating in low effective stress regimes and with overall low stress drops, as in the case of slow slip events documented in the Cascadia subduction zone, can be explained using a simple fracture propagation model with a set of parameter values consistent with laboratory data. The relative stability of ruptures governed by slip weakening and dilatancy is linked to typically slow rupture growth rates, of the order of weeks, which is primarily controlled by the long hydraulic diffusion timescales inferred for slow slip events.
Acknowledgments Discussions with Frans Aben, Dmitry Garagash, Jessica Hawthorne, Jim Rice and Rob Viesca helped shape this work. Useful suggestions and comments from Huihui Weng and an anonymous reviewer contributed to improve the manuscript. Support from the UK Natural Environment Research Council (grants NE/K009656/1 and NE/S000852/1) and from the European Research Council under the European Union's Horizon 2020 research and innovation programme (project RockDEaF, grant agreement #804685), is gratefully acknowledged. The results of this paper can be reproduced by direct implementation of the analytical formulae and numerical methods described in the main text and appendices; no new data were generated in this work.

A Elastodynamic simulations
The elastodynamic equilibrium equation for 2D, mode III ruptures is where φ(x, t) corresponds to static and dynamic contributions to stress associated with a given slip rate history on the fault. The stress in (34) must be less than (if slip rate is zero) or equal to the fault strength (if slip rate is nonzero), given by Equation (2), together with governing equations (3), (9), and (4). Numerical solutions for stress and slip rate are determined by the method outlined in Brantut et al. (2019), and only the general principle is recalled here. Time and space are discretised with constant steps. At each time step, the previous slip rate history is used to compute the stress functional φ; the spatial convolution involved in this computation is performed in the Fourier domain. Then, the strength is evaluated by computing the pore pressure distribution along the fault, approximating the convolution integral in (9) by the mid-point rule. An estimate of slip rate at the next time step is then found by equating strength with stress. The procedure is then repeated, using the average between this estimate of slip rate and that from the previous time step in the computation of the stress transfer functional and pore pressure, to arrive at the final estimate of slip rate at the next time step. The evolution of friction coefficient and porosity in space and time is tracked by updating slip at each time and space step.
Space was discretised in 4096 steps of size ∆x/L w = 1/10, where we recall that L w = µδ w /τ r , and time was discretised in 8192 steps of size ∆t = (1/2)∆x/c s . In the simulation illustrated in Figure 6, a uniform background stress equal to τ b = 1.5τ r was used, and rupture was nucleated by setting τ b = 1.01(f p /f r )τ r in a region of size 6L w around x = 0.

C Thermal pressurisation and dilatancy
Here some closed-form solutions are presented for the thermal pressurisation and dilatancy problem, with the assumption that the coefficient of friction remains constant. In that case, the solutions are sufficiently simple to allow clear physical interpretations.

C.1 Adiabatic, undrained solution
Neglecting fluid and heat diffusion, the governing equations for pore pressure and temperature simplify into two ordinary differential equations for effective stress and temperature: Straightforward integration leads to the following solution: σ (δ) = σ 0 e −δ/δc + ∆Φ max /β * δ D /δ c e −δ/δ D − e −δ/δc (43) for effective stress and for temperature. It is readily observed that the temperature asymptotically increases by a quantity (σ 0 + ∆Φ max /β * )/Λ at large slip, so that dilatancy has the effect of resetting the initial effective stress, as reported by Rice (2006).

C.2 Slip on a plane, constant slip rate solution
For Gaussian shear strain rate profiles and dilatancy proportional to strain rate, the pore pressure evolution is given in integral form by where N is the rate of dilatant fault opening, given by N/h = d∆Φ/dt. For vanishingly small h compared to both thermal and hydraulic boundary layer widths, the expression simplifies to Further assuming constant friction and slip rate, we arrive at the following integral form for shear stress: where L * = 4( √ α hy + √ α th ) 2 (ρc/f r Λ) 2 /V and D = 4α hy /f 2 r V .
The specific form (4) for dilatancy can be rewritten in terms of fault zone opening as where ∆h is the maximum opening, equivalent to ∆Φ max × h. The solution for shear stress can be found in the Laplace domain τ (s) = τ 0 L * /s whereτ denotes the Laplace transform of shear stress and s is the transformed independent variable. By standard inversion technique we find the following solution in terms of slip: Reintroducing ∆Φ max = ∆h/h, and denoting the shear stress can be rewritten with more easily recognisable quantities as In the slip on a plane, constant slip rate assumption, the temperature evolution is given by Using again Laplace transforms and after a lengthy but straightforward calculation, we arrive at The main interest of closed form expressions (51) and (53) is that they show that (1) shear stress approaches zero and (2) temperature approaches σ 0 (1 + α hy /α th )/Λ at large slip. The latter asymptote is the same as that without dilatancy (Rice, 2006). The effect of dilatancy is thus not simply to reset the initial pore pressure: since dilatancy only occurs within the fault, pore pressure will spontaneously reequilibrate with the surrounding medium, so that at large time the effect of dilatancy becomes insignificant.