Frictional weakening leads to unconventional singularities during dynamic rupture propagation

F. Paglialunga$^a$, F. Passelègue$^b$, M. Lebihain$^c$, M. Violay$^a$

$^a$LEMR, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland

$^b$Géoazur, ESEILA, Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Sophia Antipolis, France., Nice, France

$^c$Navier, Ecole des Ponts, Univ Gustave Eiffel, CNRS, Marne-la-Vallée, France

**Abstract**

Earthquakes i.e. frictional ruptures, are commonly described by singular solutions of shear crack motions. These solutions assume a square root singularity order around the rupture tip and a constant shear stress value behind it, implying scale-independent edge-localized energy. However, recent observations of large-scale thermal weakening accompanied by decreasing shear stress potentially affecting the singularity order can challenge this assumption. In this study, we replicate earthquakes in a laboratory setting by conducting stick-slip experiments on PMMA samples under normal stress ranging from 1 to 4 MPa. Strain gauges rosettes, located near the frictional interface, are used to analyze each rupture event, enabling the investigation of shear stress evolution, slip velocity, and material displacement as a function of distance from the rupture tip. Our analysis of the rupture dynamics provides compelling experimental evidence of frictional rupture driven by enhanced thermal weakening. The observed rupture fronts exhibit unconventional singularity orders and display slip-dependent breakdown work (on-fault dissipated energy). Moreover, these findings elucidate the challenges associated
with a priori estimating the energy budget controlling the velocity and final extent of a seismic rupture, when thermal weakening is activated during seismic slip.

## 1. Introduction

Frictional rupture phenomena, including natural earthquakes, are often described by singular solutions of shear crack motions (Freund, 1979; Palmer and Rice, 1973; Rice, 1980). For such cracks, the stress field at the rupture tip is described by a square root singularity ($\xi = -0.5$), constant residual stress is expected far behind the rupture tip, and the energy balance condition equates the energy release rate $G$ (i.e. rupture growth driving force) to a constant value of fracture energy $G_c$ (i.e. resistance to rupture growth). This was confirmed by experimental and numerical observations, where the onset of frictional sliding, the evolution of the rupture speed, and the rupture length were predicted by Linear Elastic Fracture Mechanics (LEFM) (Bayart et al., 2016; Kammer et al., 2015; Kammer and McLaskey, 2019; Svetlizky and Fineberg, 2014; Xu et al., 2019), and suggesting that the fracture energy controlling the dynamics of the rupture tip might be an interface property. Such an analysis often relies on the hypothesis of negligible frictional weakening far away from the rupture tip (i.e. outside of the cohesive zone).

However, it is widely recognized that fault shear stress is likely to evolve during seismic slip due to (i) velocity and slip dependencies (Marone, 1998), (ii) activation of thermal weakening processes (Di Toro et al., 2011; Hirose and Shimamoto, 2005; Rice, 2006), (iii) dilatancy inducing fluid pressure changes (Brantut, 2020; Rice and Rudnicki, 1979; Segall et al., 2010).
These changes in the residual stress behind the rupture tip could induce a slip dependency of the apparent fracture energy (nowadays more commonly called breakdown work (Tinti et al., 2005)) estimated for natural earthquakes (Abercrombie and Rice, 2005; Lambert and Lapusta, 2020), in contrast to the LEFM definition. The breakdown work $W_{bd}$ is a quantity commonly used to study the energy balance of earthquakes and is defined as an energy term including all on-fault dissipative processes $W_{bd} = \int_{0}^{D} \tau(D') - \tau_{min}dD'$, with $\tau$ the shear stress acting on the fault, $\tau_{min}$ the minimum shear stress reached on-fault, and $D$ the fault slip. It can be observed that, by definition, $W_{bd}$ is a slip-dependent quantity. It is therefore important to be aware of how possible stress weakening may affect rupture dynamics and the energy release that controls it. Such weakening phenomenon has been observed in high-velocity friction experiments (for a comprehensive review, see (Di Toro et al., 2011)). Nonetheless, friction experiments seldom exhibit elastodynamic rupture propagation due to the imposed fault slip by the experimental procedure (Chen et al., 2021). Consequently, the experimental documentation of the coupling between fault weakening and elastodynamic rupture response remains limited (Rosakis et al., 2020; Paglialunga et al., 2022).

In these regards, our recent work highlighted that a long-tailed weakening can emerge after a first rapid weakening during frictional rupture experiments (Paglialunga et al., 2022), resulting in a slip-dependent breakdown work. Despite this observation, the rupture dynamics, analyzed through LEFM, showed to be controlled by a constant fracture energy $G_c$, in agreement with previous studies (Bayart et al., 2016; Kammer et al., 2015; Kammer and
McLaskey, 2019; Svetlizky and Fineberg, 2014; Xu et al., 2019). However, analyzing such frictional ruptures in the framework of LEFM relies on the assumption of constant residual stress behind the rupture tip. The observed long-tailed weakening could call into question this assumption and limit the framework's applicability to fully describe frictional ruptures, explaining the observed mismatch between $G_c$ and $W_{bd}$ (Paglialunga et al., 2022).

Moreover, theoretical studies have shown that continuous stress weakening can modify the singularity order controlling the stress and displacement fields around the rupture tip, deviating from the square-root singularity commonly adopted in LEFM, and leading to an unconventional singularity order ($\xi \neq -0.5$) (Garagash et al., 2011; Viesca and Garagash, 2015; Brantut and Viesca, 2017; Brener and Bouchbinder, 2021b). In particular, when frictional ruptures are described by $\xi \neq -0.5$, the stress $\sigma_{ij}$ and displacement $u_i$ fields obey respectively the following scaling relationships (Brener and Bouchbinder, 2021a): $\sigma_{ij} \propto K^{(\xi)}r^{\xi}$ and $u_i \propto K^{(\xi)}r^{(\xi+1)}/\mu$, with $K^{(\xi)}$ the $\xi$- generalized dynamic stress intensity factor, $r = x - x_{tip}$ the distance from the rupture tip of a point of observation located on the fault plane, and $\mu$ the dynamic shear modulus. These lead to the following relation: $W_{bd} \propto [K^{(\xi)}]^2 r^{(1+2\xi)/\mu}$, valid for $r > x_c$, with $x_c$ the cohesive zone size (Eq. (5) from (Brener and Bouchbinder, 2021a)). From this relation, it can be easily noticed that for $\xi = -0.5$, the $W_{bd}$ dependence on $r$ completely vanishes, making the breakdown work independent of the distance from the rupture tip. This does not happen when $\xi \neq -0.5$, for which $W_{bd}$ has a direct dependence on $r$.

So far, the occurrence of such unconventional singularities during fric-
tional ruptures has not been measured at the laboratory scale. In this pa-
per, we present new data analyzed in a recently-derived theoretical frame-
work, demonstrating the first experimental evidence of strain and stress per-
turbation caused by unconventional singularities associated with velocity-
dependent frictional weakening. These experimental findings are supported
by theoretical explanations about the emergence of unconventional singular
fields during dynamic rupture.

## 2. Methods

We performed stick-slip experiments in a biaxial apparatus working in a
2D single shear configuration under an applied normal stress ranging from
1 to 4 MPa (Figure 1a). The experimental setup is the same one used and
described in (Paglialunga et al., 2022). The tested samples consist of two
polymethylmethacrylate (PMMA) blocks of dimensions (20x10x3) cm (top
block) and (50x10x3) cm (bottom block), generating, once put into contact,
an artificial fault of (20x3) cm. The external loading is imposed using two
hydraulic pumps. The normal load is applied to the top block and kept
constant while the shear load is manually increased and applied to the bottom
block inducing, once reached the fault strength, stick-slip events. Strain
gages rosettes (oriented along 45°, 90°, 135° to the fault plane), located 1 mm
away from the frictional interface, were used to compute the local strain and
stress tensors. Note that at such distance from the fault plane, the dynamic
propagation of the rupture will induce a stress perturbation (rise of stress
ahead of the rupture tip, followed by a stress drop behind it). Such a stress
perturbation will be maximum on the fault plane (y = 0+) and decrease
away from it ($y > 0$), as predicted by fracture mechanics. This feature is not observable through macroscopic measurements of strain or stress, hence the need for local ones. The strain tensor rotation was obtained through the conversion of $\varepsilon_1, \varepsilon_2, \varepsilon_3$ into $\varepsilon_{xx}, \varepsilon_{xy}, \varepsilon_{yy}$ following:

$$ \varepsilon_{xy} = \frac{\varepsilon_3 - \varepsilon_2}{2}, \qquad (1) $$

$$ \varepsilon_{yy} = \varepsilon_1, \qquad (2) $$

$$ \varepsilon_{xx} = \varepsilon_3 + \varepsilon_2 - \varepsilon_1 \qquad (3) $$

The local strain temporal evolution shows clear perturbations concurrent with stick-slips (Figure 1b). By zooming-in in time, details of the instability can be caught (Figure 1c), showing a first (main) rupture front, followed by a series of secondary fronts probably caused by rupture reflections at the fault edges. To study the rupture dynamics, only the main front was considered, selecting a time window around the first strain perturbation (Figure 1d-e). Together with the vicinity of the measurement location to the fault, the measuring interval is of paramount importance for properly addressing the elastodynamics of the studied events. From here the need for a sufficiently high recording frequency (2 MHz in this case) and a well-defined and short time interval excluding possible subsequent events. Note that the following analysis and discussions will exclusively focus on the dynamics of the main rupture front for each stick-slip event. All the experimental curves that will be shown will refer to a defined time window, systematically smaller than the expected propagation time along the fault interface (for example, the rupture showed in Figure 1e is described by a temporal window of ~45 µs). The rupture propagation velocity $C_f$ was estimated by computing the ratio between
<figure><img src="image_1.png" /><figcaption>Figure 1: a. Experimental setup - Direct shear biaxial apparatus with PMMA samples generating an artificial fault. Strain gauges rosettes are located along the fault at a distance of 1-1.5 mm from the fault plane. b. Temporal evolution of vertical strain (obtained through high-frequency strain gauges acquisition system) at the three different locations along the fault. When the fault experiences instability, the shear rupture propagates along the interface and causes a strain perturbation concurrent with the passage of the front (indicated by the blue arrows). Yellow shaded areas indicate the time window selection shown in the following panel. c. Zoom-in of (b). d. Zoom-in of (c). e. Vertical strain temporal evolution for the central location. Please note that the y-axis and x-axis limits change for each panel. t0 refers to the arrival time of the rupture for the central strain gauge.</figcaption></figure>
the distance among the strain gauge locations and the rupture front travel time from one location to the other. For each event, the particle velocity was then computed through the strain component parallel to the slip direction as $\dot{u}_x = -C_f \varepsilon_{xx}$. This estimate has been shown to be comparable to distinct measurements of slip motions associated with the propagation of the seismic rupture in previous experimental studies (Svetlizky and Fineberg, 2014; Paglialunga et al., 2022). The fault slip velocity was considered equal to twice the particle velocity measured through the strain gauges ($V = 2\dot{u}_x$), assuming an anti-symmetrical distribution of slip and slip rate. This assumption seems legitimate given that the two samples have comparable dimensions, the same width, and are made of the same material. Integrating $V$ during the propagation time, local horizontal displacement $u_x$ could be estimated as well. The slip displacement $D$ of the fault is computed as twice (refer to the assumption described just above) the horizontal displacement ($D = 2u_x$), assuming the material displacement measured through the strain gauge 1 mm away from the fault is comparable to the one occurring on-fault. The acquired temporal evolutions of strain were converted into spatial evolutions following the methods used in Svetlizky and Fineberg (2014), and by assuming that the rupture velocity obtained was locally constant during rupture propagation. Moreover, by assuming plane stress conditions, the stress tensor was computed through the elastic properties of PMMA. Given the viscoelastic nature of PMMA, dynamic values of its elastic modulus were used in this study, $E = 5.7$ GPa, and a Poisson ratio $\nu = 0.33$, both obtained through bench measurements of seismic velocities ($C_P = 2705$ m/s and $C_S = 1340$ m/s).
## 3. Results

Each rupture event was studied through the evolution of shear stress, slip velocity, and material displacement as a function of the distance from the rupture tip (Fig. 2). In such spontaneous ruptures, the slip (or slip rate) history results from the coupling of elastodynamics and the constitutive friction law, rather than being prescribed. In all the studied events, local shear stress evolution exhibited an increase ahead of the rupture tip followed by a first significant decrease within the first micrometers of slip and a second mild one within larger distances (Fig. 2a) as recently observed (Paglialunga et al., 2022). A rapid increase of slip velocity was observed concurrent with the passage of the rupture front, followed by a slow decay occurring with distance from the rupture tip. The peak slip velocity $V_{max}$ showed a clear dependence with estimated rupture speed, with ~ 0.08 m/s for $C_f$ ≈ 220 m/s up to ~ 0.8 m/s for $C_f$ ≈ 840 m/s (Fig. 2b). The evolution of horizontal displacement $u_x$ presented values close to 0µm ahead of the rupture tip (values slightly deviate from 0 due to off-fault measurement) and a sharp increase behind it (Fig. 2c), with final displacements ranging between 3.9 and 28 µm. Subsequently, the fault strength weakening was analyzed through the evolution of the local shear stress $\tau$ with the fault's slip displacement $D$. The fault's weakening presents a sharp decrease of shear stress occurring within the first microns of slip, followed by a milder decrease occurring within a larger amount of slip (Fig. 3a). The breakdown work evolution was computed as

$$W_{bd} = \int_{D((x-x_{tip})=0)}^{D} (\tau - \tau(D)) dD \quad (4)$$
<figure><img src="image_2.png" /><figcaption>Figure 2: Elastic fields around the rupture tip. Evolution of a. shear stress computed from the measured shear strain $\epsilon_{xy}$, b. slip velocity computed from the measured horizontal strain $\epsilon_{xx}$, and c. material displacement computed from the estimated slip velocity for several events presenting different $C_f$ (colorbar).</figcaption></figure>

where $D((x-x_{\text{tip}}) = 0)$ is the displacement at the passage of the rupture tip. Since no slip is expected to occur ahead of the rupture tip on the fault plane ($C_f = 0$ when $(x-x_{\text{tip}}) > 0$), the breakdown work evolution was computed only from slip occurring after the passage of the rupture tip ($x - x_{\text{tip}}) = 0$, neglecting fictitious contributions due to elastic strain of the bulk at the measurement location. The evolution of $W_{\text{bd}}$ showed a first increase with slip described by a slope close to 1 : 2 and a subsequent increase described by a slope of ~1 : 0.6(±0.1) (Fig. 3b), suggesting the existence of anomalous singularities ($\xi \neq -0.5$). Given our interest in studying the second branch of this evolution, the slope of it was measured. This was done by fitting the evolution of $W_{\text{bd}}$ with $D$ for $D > D_c$ with a first-degree polynomial in log-log scale. Then, $\xi$ was derived from estimates of the power law exponent
<figure><img src="image_3.png" /><figcaption>Figure 3: Slip-dependent breakdown work and the emergence of unconventional singularities. a. Evolution of ($\tau - \tau_{res}$) with D defining the fault's weakening for different events. The integration of these curves leads to the evolution of $W_{bd}$ with D for different $C_f$ (b). For each curve, the slope of the second branch of the evolution corresponds to a given value of singularity order. c. Evolution of estimated $\xi$ values with peak slip velocity $V_{max}$.</figcaption></figure>

through:

$$W_{bd}(D) = G_c \left( \frac{D}{D_c} \right)^{\left( \frac{1+2\xi}{1+\xi} \right)}, \quad (5)$$

for large slip displacement $D \gtrsim D_c$. The values were found to range between -0.4 and -0.2 (Fig. 3c). Another possible route proposed by Brener and Bouchbinder (2021a) is to estimate the deviation $\Delta\xi = \xi + 1/2$ of $xi$ from the LEFM value $-1/2$ from the dependence of the breakdown energy in the distance $r = x - x_{tip}$ from the crack tip. However, this requires prior knowledge of the fracture energy $G_c$ and the process zone size $x_c$, which are difficult to constrain from experimental measurements.
## 4. Theoretical modeling of the kinematic fields around the rupture tip for unconventional singularity order

While the first increase of breakdown work with slip can be explained by a slip-weakening behavior of the fault, the subsequent increase (power law of 1:0.6) is unexpected from the conventional theory of LEFM. If such a continuous weakening stage controlled the dynamics of the rupture, stress fields with a scaling $\sigma \propto r^\xi$ should be observed behind the rupture tip, as expected from theoretical studies (Brantut and Viesca, 2017; Brener and Bouchbinder, 2021b; Garagash et al., 2011; Viesca and Garagash, 2015), with the singularity order $\xi$ different from the square root singularity. To further investigate the dynamics of rupture, the temporal evolution of the strain perturbations generated by the passage of the rupture front ($\Delta\varepsilon_{xy}$, $\Delta\varepsilon_{xx}$) was compared to the theoretical predictions obtained considering both a square root singularity (LEFM) and an unconventional singularity (Brener and Bouchbinder, 2021a).

For the LEFM theoretical prediction, the stress field perturbation around the rupture tip takes the following general form (for a detailed description please refer to Freund (1998); Anderson (2017)):

$$ \Delta\sigma_{ij}(r, \theta) = \frac{K_{II}}{\sqrt{2\pi r}} \Sigma_{ij}^{II}(\theta, C_f) \qquad (6) $$

where $K_{II}$ the stress intensity factor, and $\Sigma_{ij}^{II}(\theta, C_f)$ the angular variation function. Coordinates are expressed in the polar system with $(r, \theta)$ respectively the distance from the crack tip and the angle to the crack's plane.

In the unconventional theory framework, the stress fields were derived from the elastodynamic equations assuming a steady-state rupture velocity.
The equations obtained present the following form:

$$ \sigma_{xx}(r, \theta) = \frac{4(\xi + 1)K_{\text{II}}^{(\xi)}}{\sqrt{2\pi}R(C_f)} [\alpha_s(1 + 2\alpha_d^2 - \alpha_s^2)r_d^\xi \sin(\xi\theta_d) - \alpha_s(1 + \alpha_s^2)r_s^\xi \sin(\xi\theta_s)], \quad (7) $$

$$ \tau(r, \theta) = \sigma_{xy}(r, \theta) = \tau_{\text{min}} + \frac{2(\xi + 1)K_{\text{II}}^{(\xi)}}{\sqrt{2\pi}R(C_f)} [4\alpha_s\alpha_d r_d^\xi \cos(\xi\theta_d) - (1+\alpha_s^2)^2 r_s^\xi \cos(\xi\theta_s)], \quad (8) $$

$$ \sigma_{yy}(r, \theta) = -\frac{4(\xi + 1)\alpha_s(1 + \alpha_s^2)K_{\text{II}}^{(\xi)}}{\sqrt{2\pi}R(C_f)} [r_d^\xi \sin(\xi\theta_d) - r_s^\xi \sin(\xi\theta_s)], \quad (9) $$

where $K_{\text{II}}^{(\xi)} = \lim_{r \to 0} \left( \frac{(2\sqrt{2\pi})}{(\xi+1)} r^{-\xi} \tau(r, 0^{+-}) \right)$ is the $\xi$-generalized stress intensity factor, $\alpha_d = 1 - \left(\frac{C_f}{C_d}\right)^2$ and $\alpha_s = 1 - \left(\frac{C_f}{C_s}\right)^2$, where $(C_d, C_s)$ are respectively the P-wave and S-wave velocity, and $R(C_f) = 4\alpha_d\alpha_s - (1+\alpha_s^2)^2$ is the Rayleigh function. $(r, \theta)$ are corrected for the distortion induced by the dynamic rupture velocity $C_f$, becoming $\theta_d = \arctan(\alpha_d \tan(\theta))$, $\theta_s = \arctan(\alpha_s \tan(\theta))$ and $r_d = r\sqrt{1 - \left(\frac{C_f \sin(\theta)}{C_d}\right)^2}$, $r_s = r\sqrt{1 - \left(\frac{C_f \sin(\theta)}{C_s}\right)^2}$.

The displacement field related to the unconventional rupture phenomenon can be predicted by (Brener and Bouchbinder, 2021a):

$$ u_x(r, \theta) = \frac{2K_{\text{II}}^{(\xi)}}{\mu\sqrt{2\pi}} \left[ 2\alpha_s r_d^{\xi+1} \sin[(\xi+1)\theta_d] - \alpha_s(1+\alpha_s^2)r_s^{\xi+1} \sin[(\xi+1)\theta_s] \right] \quad (10) $$

$$ u_y(r, \theta) = \frac{2K_{\text{II}}^{(\xi)}}{\mu\sqrt{2\pi}} \left[ 2\alpha_s\alpha_d r_d^{\xi+1} \cos[(\xi+1)\theta_d] - (1+\alpha_s^2)r_s^{\xi+1} \cos[(\xi+1)\theta_s] \right]. \quad (11) $$

The values of $\xi$ used to describe the experimental curves were obtained through the measured evolution of $W_{\text{bd}}$ with $D$ as discussed earlier (Brener and Bouchbinder, 2021a). The stress intensity factor was computed as (Eq. (5)) from (Brener and Bouchbinder, 2021a)): $K_{\text{II}}^{(\xi)} = \frac{EW_{\text{bd}}(D_{\text{fin}})}{(1-\nu^2)f_{\text{II}}(C_f)r_{\text{fin}}^{(1+2\xi)}}$, with $E, \nu$ respectively the elastic modulus and Poisson's ratio, and $f_{\text{II}}(C_f) = \frac{\alpha_s C_f^2}{(1-\nu)R(C_f) C_s^2}$
the universal function of rupture velocity.

## 5. Description of strain perturbations with theoretical predictions

We now compare the theoretical predictions to experimental strain and displacement evolution of two different frictional ruptures presenting values of $\xi = -0.32, -0.27$, and final values of $W_{bd}$ of 9.5 and 11 J/m$^2$, respectively (Fig. 3b). The values of singularity order were obtained as described in the Results section, by computing the slope of the second branch of the $W_{bd}$ vs $D$ curve and using Eq. 5.

The above-mentioned comparison is presented in Fig. 4. Note that for both models, i.e. LEFM and unconventional theory, the predictions of strain fail ahead of the rupture tip. This is explained by the fact that the two models assume a dynamic rupture driven along an infinite fault by a shear stress equal to the residual stress. As such, they overlook any finite-size effects emerging from the finiteness of the specimen size and the distance to the applied boundary conditions. Moreover, please note that the measurement location was chosen to be the closest possible to the fault plane (strain gauges at $\sim 1$ mm), to capture stress and displacement evolution close to the ones occurring on-fault. However, this choice implies the likelihood of performing measurements within the cohesive zone, expected to be for PMMA around 2-5 mm. The cohesive zone (indicated in Fig. 4a-d with the shaded grey area) was excluded when performing the LEFM fits, given that this model assumes conditions of small-scale yielding (dissipation zone small with respect to the other length scales).

The experimental data were compared with the predictions of LEFM
(ξ = -0.5) inverting G_c from the best possible fit. The inversion and the minimization algorithm employed to obtain the best solution of G_c use simultaneously two strain components (Δε<sub>xx</sub>, Δε<sub>xy</sub>) following the method described in previous studies (Svetlizky and Fineberg, 2014) (Fig. 4). Δε<sub>xx</sub> and Δε<sub>xy</sub> are obtained by subtracting the initial strain from ε<sub>xx</sub> and the residual strain from ε<sub>xy</sub>.

The best fits output values of G_c slightly different from the values of W_bd estimated through the integration of the slip stress curves. The LEFM predictions do not deviate excessively from the experimental curves for either event, showing an acceptable but not accurate description of the strain perturbations for Δε<sub>xx</sub> and Δε<sub>xy</sub>(Fig. 4a-b). A stronger deviation is observed for Δε<sub>xy</sub>, particularly in the case of ξ = -0.27, independently of the distance from the rupture tip (Fig. 4b). In the second stage, predictions accounting for the unconventional model were computed. The values of ξ and W_bd measured as described in the Results section were imposed. The unconventional model returned, for the two events, satisfactory predictions of the evolution of Δε<sub>xx</sub> and Δε<sub>xy</sub> (Fig. 4a,b). It can be noted that the greater the deviation from ξ = -0.5, the greater the disparities between LEFM and the unconventional model (Fig. 4a,b). In addition, the prediction obtained for u<sub>x</sub> (Fig. 4c,d) is close to the experimental evolution in terms of magnitude. However, while u<sub>x</sub> evolution is similar within the first microns, the experimental data deviate from the theoretical prediction far behind the rupture tip (Fig. 4c,d). The model returned reasonable predictions of u<sub>x</sub> for ξ=-0.32, and adequate ones for ξ=-0.27.

Finally, we compare the experimental data to both models' theoretical
<figure><img src="image_4.png" /><figcaption>Figure 4: Strain and displacement field described by unconventional singularity for two different events (respectively top and bottom panels). a., b. Comparison of the measured strain perturbations $\Delta\epsilon_{xx}$ and $\Delta\epsilon_{xy}$ with the theoretical predictions considering: i) the estimated unconventional singularities respectively $\xi=-0.32$ (a) and $\xi=-0.27$ (b), and $G = W_{bd}$ (in black) and ii) the LEFM conventional singularity $\xi=-0.5$ with $G = G_c$ (the best fit)(in grey). c., d. Evolution of the material displacement $u_x$ with predictions for unconventional and conventional singularity. e., f. Comparison of the experimental evolution of breakdown work with slip estimated at gauge location with theoretical predictions for unconventional theory (black solid line) and LEFM (grey solid line).</figcaption></figure>
299 predictions of the evolution of breakdown work with slip behind the crack tip.
300 Starting from the stress evolution estimates computed for both LEFM and
301 unconventional model, the breakdown work was computed following Eq. 4.
302 Please note that neither of the two models is expected to describe the exper-
303 imental evolution for small D, given that both of them rely on a small-scale
304 yielding condition. Concerning the overall evolution, LEFM predictions devi-
305 ate in both quantity and temporal evolution from the experimental data. On
306 the contrary, the unconventional model provides a good prediction, partic-
307 ularly for D > Dc, as expected from the unconventional theory (Fig. 4e-f).
308 These results highlight that while LEFM provides reasonable estimates of
309 fracture energy, the unconventional theory provides more coherent predic-
310 tions of breakdown work evolution with slip when enhanced weakening is
311 observed.

6. Flash heating as a possible weakening mechanism

313 These results provide the first complete evidence of unconventional stress
314 fields during the dynamic propagation of laboratory frictional rupture, caused
315 by continuous stress weakening behind the rupture tip. The observed un-
316 conventional singularity orders could emerge, among others, from frictional
317 weakening mechanisms such as; thermal activation (Bar-sinai et al., 2014),
318 viscous friction (Brener and Marchenko, 2002), powder lubrication (Reches
319 and Lockner, 2010), flash heating (Molinari et al., 1999; Rice, 2006; Brantut
320 and Viesca, 2017), thermal pressurization (Rice, 2006; Viesca and Garagash,
321 2015). Among these, flash heating has been shown to be activated under
322 similar experimental conditions in Homalite (Rubino et al., 2017), and thus
could be the best candidate to explain the unconventional stress fields observed in our experiments. Moreover, the high slip rate measured near-fault promotes the activation of flash heating as previously shown (Molinari et al., 1999; Rice, 2006; Goldsby and Tullis, 2011). This agrees with the clear dependence of $\xi$ values with maximum slip rate and rupture velocity observed in our events (Fig. 3c): higher $V_{max}$ are associated with $\xi$ values that deviate from the conventional value (-0.5).

Flash heating is activated when the fault slip velocity becomes higher than a critical weakening slip velocity $V_w$, causing mechanical degradation of contact asperities during their lifetime (Rice, 2006; Goldsby and Tullis, 2011). The temperature reached at the asperities was computed through $T_{asp} = T_{amb} + \frac{1}{(\rho c_p \sqrt{k \pi})} \tau_c V \sqrt{t_c}$ with $T_{amb}$ the ambient initial temperature, $\tau_c$ the stress acting on the single asperity, $t_c$ the lifetime of a contact, $\rho$ the bulk density, $c_p$ the bulk specific heat and $k$ the thermal diffusivity. Under our experimental conditions, the temperature increased with slip velocity, exceeding the material's melting temperature ($T_{asp} > T_{melting} = 160^{\circ}\text{C}$) (Fig. 5a, b), and indicating that melting of asperities probably occurred in our experiments (Rubino et al., 2017). We compared the evolution of $W_{bd}$ with $D$, normalized respectively by $G_c$ and $D_c$, with asymptotic solutions for flash heating phenomena (Brantut and Viesca, 2017).

For $D < D_c$ (small slip), the evolution of $W_{bd}$ can be described by the asymptotic solution derived for adiabatic conditions (Brantut and Viesca, 2017):

$$W_{bd} = \rho c (T_m - T_{amb}) w \sqrt{2\pi} \left( \frac{D}{V t_w^A + D} \right)^2 \quad (12)$$

where $t_w^A = \rho c (T_m - T_f) / \tau_a (\sqrt{2\pi w}) / V_w$ (time required for a layer of thickness
347 $\sqrt{2\pi w}$ to reach $T_{\text{melting}}$), $w$ is the shear zone thickness (assumed here as 348 $w = 4a$ with $a$ the asperity size), and $\tau_a$ is a normal stress dependent contact 349 shear stress at the origin of the change in temperature in the fault layer 350 (Fig. 5c). In presence of gouge along the interface, $\tau_a$ will correspond to 351 the macroscopic shear stress $\tau_0$. Along bare rock interfaces, $\tau_a = \tau_c \frac{a}{\Delta L_{\text{asp}}}$, 352 where $\Delta L_{\text{asp}}$ is the average distance between two asperities (see Annex A for 353 details). Note that this model assumes a constant sliding velocity $V$. This 354 assumption looks fairly reasonable in our case, as the first part of the stress 355 weakening ($D < D_c$) occurs in a very short time window, during which $V$ is 356 nearly constant.

357 For $D > D_c$, a second asymptotic solution considering the coupled elas- 358 todynamics and frictional motions of the propagating rupture can be used 359 (Brantut and Viesca, 2017):

$$W_{\text{bd}} = \tau_c D_w^{SP} \left( \frac{\mu V_w}{3\pi \tau_a C_f} \right)^{(1/3)} \left( \frac{D}{D_w^{SP}} \right)^{(2/3)} \quad (13)$$

360 where $D_w^{SP} = V_w \alpha (\frac{\rho c(T_w - T_f)}{\tau_a V_w})^2$ is a characteristic slip weakening distance. 361 While this asymptotic solution is expected to describe the evolution of break- 362 down work at a larger seismic slip than the one observed in our experiments, 363 this equation can still be used here because (i) heat diffusion at the scale of 364 asperities is expected to control fault weakening when $D > D_c$ and (ii) $\tau_a$ 365 increases with $\tau_0$, through the increase of $\frac{a}{\Delta L_{\text{asp}}}$ with $\sigma_n$. 366

366 Assuming our experimental estimate of $C_f$, this asymptote well describes the 367 second branch of the evolution of $W_{\text{bd}}$ with $D$ (power law with an exponent 368 of 2/3, Fig. 5c). Such scaling is also observed at large slip for thermal pres- 369 surization in drained conditions, suggesting that this exponent is related to
<figure><img src="image_5.png" /><figcaption>Figure 5: a. Evolution of local shear stress $\tau$, with slip velocity for one event. b. Temperature evolution with slip velocity at asperity scale compared with melting temperature of PMMA ($T_m$ = 160). c. Slip dependence of breakdown work (curves are normalized respectively by $G_c$ and $D_c$). $W_{bd}$ evolution exhibits two power laws with exponents of 2 and 0.6. The experimental curves are all described by the asymptotic solutions related to an adiabatic regime for small $D$ and a diffusive regime for large $D$ (Brantut and Viesca, 2017). The dotted black line shows the expected evolution of $W_{bd}$ assuming LEFM at the strain gauges position.</figcaption></figure>
diffusion mechanisms regulating the weakening of faulting during seismic slip (Brantut and Viesca, 2017; Viesca and Garagash, 2015).

Importantly, an energy dissipation $W_{bd}$ greater than the fracture energy $G_c$ was already observed in Barras et al. (2020) for sliding interfaces whose frictional behavior is described by a rate-and-state friction law. Despite this excess, the rupture dynamics were well described by a conventional LEFM analysis (with $\xi = -0.5$). This was later justified by Brener and Bouchbinder (2021b), who showed that ruptures along interfaces obeying rate-and-state friction displayed a singularity $\xi = -0.406 \simeq -0.5$, which corresponds to the lower end of our measurements. However, fault characteristics (e.g. roughness, fluid diffusivity, etc.) and external factors such as initial stress state or on-fault temperature can alter the friction law that controls interface slip (i.e. flash heating, thermal pressurization, and others) and change the singularity observed near the rupture accordingly. In the case of flash heating, the observed evolution of breakdown work with slip generates, for example, a singularity order $\xi = -0.25$ (Brantut and Viesca, 2017), which corresponds to the higher-end exponents of Fig. 3. In our experiments, continuous values of exponents $\xi$ have been measured between $\xi = -0.42$ (rate-and-state) and $\xi = -0.22$ (flash heating). This can be caused by the presence of a population of contact asperities, each of which have a different size, experience a different normal and shear stress, and reach thus a different value of temperature and slip velocity during rupture (implying that not necessarily all contact asperities experience flash heating) resulting, on average, in a smooth transition from rate and state frictional contact for the lower slip velocities (nearly conventional, $\xi = -0.5$) to flash heating for larger slip
395 velocities (unconventional, $\xi = -0.25$).

396 7. Implications and conclusions

397 These experimental results show that the continuous weakening activated
398 along the fault can modify the singularity order governing displacement and
399 stress fields around the rupture tip, inducing a slip and scale-dependent
400 breakdown work, rather than a constant one. Moreover, this work high-
401 lights from an experimental point of view that frictional rupture analysis in
402 the linear elastic fracture framework might not always be sufficient when fric-
403 tional weakening mechanisms occur away from the rupture tip. Importantly,
404 as long as the residual stress does not reach a steady-state value far from
405 the rupture tip, as happens for thermal weakening processes, the singular
406 fields will hardly recover the conventional square-root singularity, indepen-
407 dently of the rupture size. One could nonetheless assess the dynamics of
408 such earthquakes, building on a Griffith criterion adapted to unconventional
409 singularities (see Eq. (7) (Brener and Bouchbinder, 2021a)). However, this
410 would involve both the fracture energy and the cohesive zone size that of-
411 ten depends on the structural problem (loading conditions, fault geometry).
412 Furthermore, the activation of thermal mechanisms depends not only on the
413 rupture characteristics such as crack velocity but also on ambient conditions
414 (such as initial temperature) and possibly slip history controlling asperity
415 roughness and strength. As a result, both rupture dynamics and fault weak-
416 ening are expected to be governed by fault geometry and rheology and may
417 vary depending on the natural environment.

418 Our new results highlight the difficulty in a priori estimating the relevant
parameters governing the dynamics of the seismic rupture, expected to con-
trol the final rupture length (earthquake size). One may legitimately wonder
whether theoretical models will be able to capture these complex behaviors,
or whether numerical simulations, as proposed in recent studies, will be re-
quired instead (Lambert and Lapusta, 2020).

However, together with the recent development of the unconventional singu-
larity theory (Brener and Bouchbinder, 2021a), our results open the door for
a better understanding of the rupture dynamics and energy budget of natural
earthquakes, through the possible evaluation of the equations of motions for
unconventional rupture phenomena.
## 429 Appendix A.

For the estimate of $\Delta L_{asp}$, a simplified description of the interface roughness is used, considering only one population of asperities of typical size a and height h, separated by an average distance $\Delta L_{asp}$. The number of asperities was computed considering the following relationship $\frac{A_r}{A_n} = \frac{G_c}{G_{PMMA}}$ (values of $G_{PMMA}$ coming from Marshall et al. (1974)), which lead to $N_{2D} = \frac{G_c}{G_{PMMA}} \frac{A_n}{\pi D_{asp}^2/4}$.

Assuming an equidistant spacing between the asperities in both directions, the total number of asperities can be written as $N_{2D} = N_x N_y$ with $N_x$ and $N_y$ respectively the number of rows and columns of asperities located in the x and y directions. The latter numbers are related to the interface dimensions through $\frac{N_x}{N_y} = \frac{L_f}{W_f}$, with $L_f$ and $W_f$ respectively the length and width of the interface. Considering this as a 1-D problem, the number of asperities along the interface in the slip direction reads $N_{1D} = \sqrt{N_{2D} \frac{L_f}{W_f}}$. The distance between two asperities could then be estimated as $\Delta L_{asp} = \frac{L_f - N_{1D} D_{asp}}{N_{1D} + 1}$. The contact stress at the origin of the change in temperature of asperities during the seismic slip can be expressed as $\tau_{eff} = \tau_c \frac{a}{\Delta L_{asp}}$.

## 446 Appendix B.

Acknowledgments.

F.P. and M.V. acknowledge support by the European Research Council Starting Grant project 757290-BEFINE.
450 References

451 Abercrombie, R.E., Rice, J.R., 2005. Can observations of earthquake scaling
constrain slip weakening? Geophysical Journal International 162, 406–424.
doi:10.1111/j.1365-246X.2005.02579.x.

454 Anderson, T.L., 2017. Fracture mechanics: fundamentals and applications.
CRC press.

456 Bar-sinai, Y., Spatschek, R., Brener, E.a., Bouchbinder, E., 2014. On the
velocity-strengthening behavior of dry friction. Journal of Geophysical
Research 119, 1738–1748. doi:10.1002/2013JB010586.Abstract.

459 Barras, F., Aldam, M., Roch, T., Brener, E.A., Bouchbinder, E., Molinari, J.F., 2020. The Emergence of Crack-like Behavior of Frictional Rupture: Edge Singularity and Energy Balance. Earth and Planetary Science Letters 531, 115978. doi:10.1016/j.epsl.2019.115978.

463 Bayart, E., Svetlizky, I., Fineberg, J., 2016. Fracture mechanics determine the lengths of interface ruptures that mediate frictional motion. Nature Physics 12, 166–170. doi:10.1038/nphys3539.

466 Brantut, N., 2020. Dilatancy-induced fluid pressure drop during dynamic rupture: Direct experimental evidence and consequences for earthquake dynamics. Earth and Planetary Science Letters 538, 116179. URL: https://doi.org/10.1016/j.epsl.2020.116179, doi:10.1016/j.epsl.2020.116179, arXiv:1904.10906.

471 Brantut, N., Viesca, R.C., 2017. The fracture energy of ruptures
driven by flash heating. Geophysical Research Letters 44, 6718-6725.
doi:10.1002/2017GL074110.

Brener, E., Bouchbinder, E., 2021a. Theory of unconventional singularities
of frictional shear cracks. Journal of the Mechanics and Physics of Solids
153, 104466. doi:10.1016/j.jmps.2021.104466.

Brener, E., Bouchbinder, E., 2021b. Unconventional singularities and energy balance in frictional rupture. Nature Communications 12, 2585.
doi:10.1038/s41467-021-22806-9.

Brener, E.A., Marchenko, V.I., 2002. Frictional shear cracks. JETP Letters
76, 211-214. doi:10.1134/1.1517386, arXiv:0204046.

Chen, X., Chitta, S.S., Zu, X., Reches, Z., 2021. Dynamic fault weakening
during earthquakes: Rupture or friction? Earth and Planetary Science
Letters 575, 117165.

Di Toro, G., Han, R., Hirose, T., De Paola, N., Nielsen, S., Mi-
zoguchi, K., Ferri, F., Cocco, M., Shimamoto, T., 2011. Fault
lubrication during earthquakes. Nature 471, 494-499. URL:
http://dx.doi.org/10.1038/nature09838, doi:10.1038/nature09838.

Freund, L.B., 1979. The mechanics of dynamic shear crack propagation. Journal of Geophysical Research: Solid Earth 84, 2199-2209.
doi:10.1029/JB084iB05p02199.

Freund, L.B., 1998. Dynamic fracture mechanics. Cambridge university
press.
Garagash, D.I., Detournay, E., Adachi, J.I., 2011. Multiscale tip asymptotics in hydraulic fracture with leak-off. Journal of Fluid Mechanics 669, 260–297. doi:10.1017/S002211201000501X.

Goldsby, D., Tullis, T., 2011. Flash Heating Leads to Low Frictional Earth-quake Slip Rates. Science 334, 216–218.

Hirose, T., Shimamoto, T., 2005. Growth of molten zone as a mechanism of slip weakening of simulated faults in gabbro during frictional melting. Journal of Geophysical Research: Solid Earth 110, 1–18. doi:10.1029/2004JB003207.

Kammer, D.S., McLaskey, G.C., 2019. Fracture energy estimates from large-scale laboratory earthquakes. Earth and Planetary Science Letters 511, 36–43. doi:10.1016/j.epsl.2019.01.031.

Kammer, D.S., Radiguet, M., Ampuero, J.P., Molinari, J.F., 2015. Linear elastic fracture mechanics predicts the propagation distance of frictional slip. Tribology Letters 57. doi:10.1007/s11249-014-0451-8, arXiv:1408.4413.

Lambert, V., Lapusta, N., 2020. Rupture-dependent breakdown energy in fault models with thermo-hydro-mechanical processes. Solid Earth 11, 2283–2302. doi:10.5194/se-11-2283-2020.

Marone, C., 1998. Laboratory-derived friction laws and their application to seismic faulting. Annual Review of Earth and Planetary Sciences 26, 643–696. doi:10.1146/annurev.earth.26.1.643.
516 Marshall, G., Coutts, L.H., Williams, J., 1974. Temperature effects in the
517 fracture of pmma. Journal of Materials Science 9, 1409-1419.

518 Molinari, A., Estrin, Y., Mercier, S., 1999. Dependence of the coefficient of
519 friction on the sliding conditions in the high velocity range. Journal of
520 Tribology 121, 35-41. doi:10.1115/1.2833808.

521 Paglialunga, F., Passelègue, F., Brantut, N., Barras, F., Lebihain, M., Violay,
522 M., 2022. On the scale dependence in the dynamics of frictional rupture
523 : Constant fracture energy versus size-dependent breakdown work. Earth
524 and Planetary Science Letters 584. doi:10.1016/j.epsl.2022.117442.

525 Palmer, A.C., Rice, J.R., 1973. The growth of slip surfaces in the pro-
526 gressive failure of over-consolidated clay. Royal society 548, 527-548.
527 doi:10.1098/rspa.1973.0040.

528 Reches, Z., Lockner, D.A., 2010. Fault weakening and earthquake instability
529 by powder lubrication. Nature 467, 452-455. doi:10.1038/nature09348.

530 Rice, J.R., 1980. The mechanics of earthquake rupture. Physics of the Earth's Interior , 555-649URL: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.161.3251, doi:10.1.1.161.3251.

531 Rice, J.R., 2006. Heating and weakening of faults during earthquake slip. Journal of Geophysical Research: Solid Earth 111, 1-29. doi:10.1029/2005JB004006.

532 Rice, J.R., Rudnicki, J.W., 1979. Earthquake precursory effects due to pore
fluid stabilization of a weakening fault zone. Journal of Geophysical Research: Solid Earth 84, 2177–2193. doi:10.1029/JB084iB05p02177.

Rosakis, A.J., Rubino, V., Lapusta, N., 2020. Recent Milestones in Unraveling the Full-Field Structure of Dynamic Shear Cracks and Fault Ruptures in Real-Time: From Photoelasticity to Ultrahigh-Speed Digital Image Correlation. Journal of Applied Mechanics 87, 030801. doi:10.1115/1.4045715.

Rubino, V., Rosakis, A.J., Lapusta, N., 2017. Understanding dynamic friction through spontaneously evolving laboratory earthquakes. Nature Communications 8. doi:10.1038/ncomms15991.

Segall, P., Rubin, A.M., Bradley, A.M., Rice, J.R., 2010. Dilatant strengthening as a mechanism for slow slip events. Journal of Geophysical Research: Solid Earth 115, 1–37. doi:10.1029/2010JB007449.

Svetlizky, I., Fineberg, J., 2014. Classical shear cracks drive the onset of dry frictional motion. Nature 509, 205–208. doi:10.1038/nature13202.

Tinti, E., Spudich, P., Cocco, M., 2005. Earthquake fracture energy inferred from kinematic rupture models on extended faults. Journal of Geophysical Research: Solid Earth 110, 1–25. doi:10.1029/2005JB003644.

Viesca, R.C., Garagash, D.I., 2015. Ubiquitous weakening of faults due to thermal pressurization. Nature Geoscience 8, 875–879. doi:10.1038/ngeo2554.

Xu, S., Fukuyama, E., Yamashita, F., 2019. Robust Estimation of Rupture
560 Properties at Propagating Front of Laboratory Earthquakes. Journal of Geophysical Research 124, 766-788. doi:10.1029/2018JB016797.

561