# Laboratory Earthquakes Simulations-Typical Events, Fault Damage, and Gouge Production

Guilhem Mollon, Jerome Aubry, Alexandre Schubnel

► To cite this version:

Guilhem Mollon, Jerome Aubry, Alexandre Schubnel. Laboratory Earthquakes Simulations-Typical Events, Fault Damage, and Gouge Production. Journal of Geophysical Research: Solid Earth, 2023, 128 (2), 10.1029/2022JB025429. [hal-04661769](https://hal.science/hal-04661769v1)

HAL Id: hal-04661769

https://hal.science/hal-04661769v1

Submitted on 26 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2022JB025429

**Key Points:**

* A coupled Discrete-Continuum model of a laboratory earthquakes experiment leads to the spontaneous development of numerous seismic cycles

* The stress fields on the fault are found to be very heterogeneous. Typical events of various sizes are analyzed in details

* Damage and wear of the interface are quantified during sliding, and their rates are found to follow different laws

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

G. Mollon, guilhem.mollon@insa-lyon.fr

**Citation:**

Mollon, G., Aubry, J., & Schubnel, A. (2023). Laboratory earthquakes simulations—Typical events, fault damage, and gouge production. *Journal of Geophysical Research: Solid Earth*, 128, e2022JB025429. https://doi.org/10.1029/2022JB025429

Received 19 AUG 2022
Accepted 23 JAN 2023

**Author Contributions:**

**Conceptualization:** Guilhem Mollon, Jérôme Aubry, Alexandre Schubnel

**Data curation:** Guilhem Mollon

**Investigation:** Guilhem Mollon

**Methodology:** Guilhem Mollon

**Software:** Guilhem Mollon

**Validation:** Jérôme Aubry, Alexandre Schubnel

**Writing – original draft:** Guilhem Mollon

**Writing – review & editing:** Jérôme Aubry, Alexandre Schubnel

© 2023. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

## Laboratory Earthquakes Simulations—Typical Events, Fault Damage, and Gouge Production

Guilhem Mollon¹ ID, Jérôme Aubry² ID, and Alexandre Schubnel²

¹Université de Lyon, INSA Lyon, CNRS, LaMCoS, Villeurbanne, France, ²Laboratoire de Géologie, École Normale Supérieure/CNRS UMR 8538, PSL University, Paris, France

**Abstract** We propose a numerical model of laboratory earthquake cycle inspired by a set of experiments performed on a triaxial apparatus on sawcut Carrara marble samples. The model couples two representations of rock matter: rock is essentially represented as an elastic continuum, except in the vicinity of the sliding interface, where a discrete representation is employed. This allows to simulate in a single framework the storage and release of strain energy in the bulk of the sample and in the loading system, the damage of rock due to sliding, and the progressive production of a granular gouge layer in the interface. After independent calibration, we find that the tribosystem spontaneously evolves toward a stick-slip sliding regime, mimicking in a satisfactory way the behavior observed in the lab. The model offers insights on complex phenomena which are out of reach in experiments. This includes the variability in space and time of the fields of stress and effective friction along the fault, the progressive thickening of the damaged region of rock around the interface, and the build-up of a granular layer of gouge accommodating shear. We present in detail several typical sliding events, we illustrate the fault heterogeneity, and we analyze quantitatively the damage rate in the numerical samples. Some limitations of the model are pointed out, as well as ideas of future improvements, and several research directions are proposed in order to further explore the large numerical data set produced by these simulations.

**Plain Language Summary** Earthquakes are due to sudden sliding in faults several kilometers in the ground. A common laboratory practice is to reproduce such sliding events on dedicated lab devices. In this work, we present a novel numerical model aiming to reproduce such experiment in a computer simulation, in order to enhance our understanding of the phenomena at stake. This model is novel because it couples two different representations of the rock matter, namely a continuous and a discrete one. It therefore allows to reproduce in the same framework the bulk deformation of rock and the granular phenomena occurring at the sliding interface. The model is calibrated and leads to the spontaneous occurrence of unstable sliding, that is, of earthquakes of the same kind as those observed in the lab. We further explore into more detail some typical sliding events, and focus our attention of the interface damaging and wear during sliding. This work is likely to clarify our interpretations of sliding events in the lab.

## 1. Introduction

Laboratory earthquakes are experimental models of seismic events in the lab. In the past 60 years, they have been performed in a wide variety of experimental configurations and for various purposes. A first major class of such experiments consists in driving a controlled sliding motion between two surfaces of rock in order to reproduce realistic kinematics of a seismic event and to study the frictional response of the interface (Di Toro, 2006; Goldsby & Tullis, 2011; Hirose & Shimamoto, 2005; Reches & Lockner, 2010; Sone & Shimamoto, 2009). While most such studies focus on frictional weakening (Chang et al., 2012), some of them are also dedicated to the quantification of wear kinetics (Boneh et al., 2014) or of breakdown work (Nielsen et al., 2016). Another class of experiments tries to observe spontaneous sliding events in a preloaded interface, in order to understand the phenomena at stake during nucleation and propagation of earthquakes. This includes a large body of experimental work performed on analogous materials (Bayart et al., 2018; Latour et al., 2013; Rubino et al., 2014, 2017; Xia et al., 2004), which present the advantage of easing the failure observation and of allowing more direct quantitative measurements. But this class also includes a large number of experiments on more realistic rock samples (Aubry et al., 2020; Goebel et al., 2013; Kandula et al., 2019; Li & Zhou, 2021; Marty et al., 2019; Passelègue, Schubnel, et al., 2016; Sobolev et al., 1996), with a recent trend toward meters-size rock samples allowing a mapping of the rupture front and of its associated stress field during its progression (McLaskey, 2019; Xu et al., 2019). It is worth mentioning as well the numerous experiments performed on granular materials,
either analogous (Houdoux et al., 2021; Lherminier et al., 2019) or natural (Bolton et al., 2020; Rouet-Leduc et al., 2017), in order to study the failure properties and statistics of the layer of granular fault gouge usually found in seismic faults. At the core of these experiments is the concept of stick-slip, and the idea that the earthquake cycle is basically the same phenomenon (albeit at a much larger scale in space and time) as the high-frequency succession of adhesion and sliding episodes often observed and documented in mechanical interfaces (Feeny et al., 1998; Ferrer et al., 2010; Yoshizawa & Israelachvili, 1993).

In order to investigate the wide range of possible sliding behavior that may occur at brittle-ductile transition, it is necessary to submit rock samples to stress levels comparable to those existing at depths between 15 and 25 km, where this transition is found to occur. This corresponds to confining stresses of the order of several hundreds of MPa, while most experimental systems used to trigger and observe laboratory earthquakes are limited to a few tens of MPa at most. One notable exception is the case of triaxial cells, which can approach the appropriate stress level (Aubry et al., 2020), and can therefore put rock samples in a mechanical state close to that of the deep seismogenic zone. Measurements in such apparatus are however more challenging than in direct-shear configurations, because of the cylindrical shape of the samples and of the confining system applying the larger confining stress. Acoustic sensors and strain gauges can be deployed concurrently to detect and characterize the rupture (Marty et al., 2019), but it remains difficult to image it with a good resolution both in space and time. One solution to improve our understanding of the mechanical processes at stake in such a complex environment is to use numerical simulation. Various numerical techniques are used in the field of numerical tribology to understand the local response of a contact interface during sliding (Renouf et al., 2011; Vakis et al., 2018). Since tribological contacts nearly always fit in the tribological triplet framework (Berthier, 1990), a representation of the third body is needed. We call third body the layer of matter which fills the interface, separates the surfaces, transmits the load, and accommodates the relative displacement. The most common numerical method to represent the third body in tribological simulations is the Discrete Element Method (DEM), which consists in applying the equations of motion to each grain composing the third body, and to have these grains interact through contacts with properly chosen contact models (Cao et al., 2011; Fillot et al., 2007; Iordanoff et al., 2005; Kabir et al., 2008; Mollon, 2015, 2019; Quacquarelli, 2021; Renouf et al., 2011). The geological third body is the fault gouge, which originates from the degradation and wear of the host rock during sliding. Several authors have used DEM or FDEM (combined Finite-Discrete Element Method) to represent fault gouge, simulate its seismic shearing, and investigate its frictional response as a third body (Aharonov & Sparks, 2004; Casas et al., 2022; Cho et al., 2008; Da Cruz et al., 2005; Dorostkar et al., 2017; Gao et al., 2018; Guo & Morgan, 2004; Mair & Marone, 1999; Mair et al., 2002; Mollon et al., 2021; Morgan, 1999; Morgan & Boettcher, 1999; Zhao et al., 2012), following the same logic as in the tribological literature. Another area relevant to earthquakes mechanics that has been investigated with DEM is the fracturing of intact rock. Various approaches inspired by the Bonded-Particle-Model (Potyondy & Cundall, 2004) have been proposed in order to represent in the same numerical framework an initially elastic medium which can progressively damage and fracture (Asadi et al., 2012; Hazzard & Young, 2000; Kim & Buttlar, 2009; Park & Min, 2015; Potyondy, 2007; Wanne & Young, 2008). A recent trend in rock mechanics is to enrich the traditional circular/spherical representation of grains by considering an initial assembly of polygonal bodies (Saksala & Jabareen, 2019; Z. Wang et al., 2021). This geometry allows to consider samples without any initial porosity, and to represent more accurately the nucleation and the development of fractures and microcracks and the associated dilatancy. A similar approach was used with elastic triangular elements in Okubo et al. (2019) in order to study the off-fault damage generated by the dynamic stress field at the seismic rupture front, using the Finite-Discrete Element Methods (FDEM).

In this paper, we wish to extend these numerical methods to build a model that could reproduce numerically the seismic cycle as observed in the lab. We use the laboratory earthquake experiments of Aubry et al. (2020) as an inspiration, although we cannot yet claim that the model accurately reproduces the experimental results (for several reasons that we develop in Section 7). These experiments, that are described in more details in Section 3, consist in applying an axial strain rate to sawcut marble samples in a triaxial cell under a confining stress of tens to hundreds of MPa. In order to avoid making unnecessary assumptions concerning the friction law controlling the sliding interface, we wish the surfaces to damage along the seismic cycles and the interface to get filled with a certain layer of granular gouge controlling friction, as observed experimentally (Aubry et al., 2020). We therefore need a model where rock is able to (a) behave elastically in order to dynamically store and release strain energy along seismic cycles, (b) get damaged in the vicinity of the fault surface, (c) separate from the host rock in order to form a granular third body, that will in turn control the frictional response of the interface. This requires a
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Calibration model for a discrete model of Carrara marble, inspired by Fredrich et al. (1989). (a) full view of the simulation and boundary conditions; (b) close-in view on the upper-right corner, zooming on the displacement-controlled upper rigid plate (gray), on the continuous flexible membrane transmitting the confining pressure (green), and on the cohesive polygonal elements representing the elastic-degradable marble (yellow); (c) illustration of the normal gap between a node and a segment from two different grains; (d) cohesive Zone Model contact law (similar to Quacquarelli, 2021). Note that the mesh of the continuum part in b is only provided for graphic representation, but does not represent the nodes connectivity since meshfree shape functions are used.</figcaption></figure>

coupling between a discrete and a continuum-based representation of solid matter. In Section 2, we propose and calibrate independently a discrete model for the elastic-brittle behavior of marble. Section 3 is dedicated to the description of the reference experiment and of the discrete-continuum numerical model inspired by it. Section 4 describes the frictional response of the system and the spontaneous occurrence of seismic cycles, and Section 5 provides some insight into typical sliding events observed during simulations. Section 6 adds several quantitative results about the kinetics of damaging and wear of the surfaces in relation with the development of the seismic cycles, and Section 7 discusses these results and proposes future research directions.

## 2. A Discrete Model for Marble

Section 3 presents a Discrete-Continuum model for laboratory earthquakes, which is then more deeply analyzed in further sections. Before this presentation, though, it is necessary to design and to calibrate an appropriate model for Carrara marble. As we elaborate later, this model should be able to represent (a) the elastic response of marble at low strain levels, (b) its damaging beyond a certain yield stress, and (c) its fracturing and fragmentation under extreme stress levels. The first objective could easily be met by any classical continuum-based method, but the last two require a different approach, based on a more discrete representation of matter. In order to test and calibrate such a model, the classical experimental results from Fredrich et al. (1989) are taken as a reference. In these experiments, cylindrical samples of intact marble (diameter 15.88 mm, height 38.10 mm) were submitted to constant confining stresses $\sigma_3$ ranging from 5 to 300 MPa and to a constant compressive vertical strain rate $\epsilon_1 = 10^{-5} s^{-1}$, up to a few percents of vertical strain. During experiments, vertical stress $\sigma_1$ and volumetric strain $\epsilon_v$ were monitored. Various analysis tools were then used to investigate the micromechanics behind the measured macroscopic mechanical response of the samples.

A numerical reproduction of this experiment is represented in Figure 1a. This reproduction is at scale, but is restricted to two dimensions with plane strain kinematics. Intact rock is represented as a collection of polygonal elements, which are conformal in their initial state (meaning that porosity is negligible). This collection of polygons forms a rectangle which is surrounded by two elastic membranes (left and right, which are used to transmit the confining pressure without restricting the kinematic freedom of the sample) and two rigid plates
(above and below, which transmit the axial strain), Figure 1b. This model is introduced in the code MELODY (Mollon, 2016, 2018), which can deal both with large collections of rigid bodies with arbitrary shapes (in the DEM framework) and with highly deformable elastic bodies (in the Multibody Meshfree framework, using Moving Least Square shape functions (Belytschko et al., 1994)). The elastic deformations of the membranes follow the equations of continuum mechanics (with a Young's modulus of 1 MPa), and interact with the rigid grains through a Coulomb contact law with a friction coefficient of 0.3. The rigid grains interact through a contact law belonging to the general class of the Cohesive Zone Models (CZM, (Kim & Buttlar, 2009)). CZM are used for various applications when initially bonded objects, either rigid or deformable, are submitted to damage and eventually to separation: fracture, delamination, tribological wear, etc. In the initial stage of the simulation, all contacts between grains receive the status “intact,” along with a “Damage” parameter *D* equal to 0. Following the classical DEM procedure, at any time in the simulation, a normal and a tangential gap are computed at each contact between rigid grains. Close pairs of grains are first listed using broad proximity detection algorithms (Mollon, 2018), and gaps are computed for any node of a given body close enough to a segment of another body. More specifically, a normal gap $\delta_n$ is computed by orthogonally projecting this node on this segment (Figure 1c), and a tangential gap $\delta_t$ is obtained by integrating in time the relative tangential displacement of the two bodies. From these gaps and the current value of the bond damage, a normal and a tangential contact stress (respectively $P_n$ and $P_t$) are computed, using the CZM model described in Figure 1d. This model relies on five physical parameters: a contact stiffness *k*, an ultimate strength $P_{lim}$, a residual strength $P_{res}$, a full-damage gap $\delta_{res}$, and a contact loss gap $\delta_{loss}$. The same law is applied in the tangential and normal directions. If the damage *D* is equal to 1, the bond has the status “broken.” In that case, it is submitted to a purely cohesive contact law with a cohesion $P_{res}$. Conversely, if the damage *D* is strictly lower than 1, two cases may occur:

1. If the gap $\delta$ (either normal or tangential) is lower than its maximum value in the past history of the bond, called $\delta_{max}$ and related to the current bond damage, the bond behaves elastically.

2. If the gap $\delta$ is larger than $\delta_{max}$, the bond damage is increased, and the value of $\delta_{max}$ is updated to the current value of $\delta$.

It should be noted that if the bond is intact ($D=0$), as in the initial state of the system, the value of $\delta_{max}$ is set to $\delta_{lim} = P_{lim}/k$. As shown in Figure 1d, when the damage parameter increases, the bond stiffness and strength both decrease. The case of a broken bond ($D=1$) is reached when $\delta$ reaches the value of the parameter $\delta_{res}$, which therefore controls the amount of ductility in the bond. The contact is lost if the normal gap goes beyond the parameter $\delta_{loss}$. Any newly established contact is then initialized with a damage parameter equal to $D=1$. This means that damage and fracture are irreversible and that no healing is implemented. From the damage *D* of each bond, a grain-related damage parameter, called “Relative damage” in the remainder of this study, is computed as the mean value of the damages on all the nodes of the grain contour, weighted by the half-lengths of the two segments surrounding this node. This process allows to assign to each grain an average damage parameter (comprised between 0 for no damage, and 1 for a complete damage), which quantifies the proportion of its initially intact bonds which were damaged and eventually broken. This contact law is completed by a classical viscous damping $\gamma$ to avoid numerical instability related to the accumulation of kinetic energy in the system. The integration of all the resulting contact forces on the contour of each grain provides a resulting force and a resulting torque, which are then plugged into an explicit solver to integrate in time the equations of motion of each object. More details on the contact detection and solving in MELODY are provided in (Mollon, 2016, 2018).

A numerical calibration is first performed in order to set the grain size and the loading rate in such a way that the strength of the sample is not affected by the coarseness of the discretization and by inertial effects, respectively. Four confining stresses are applied: 5, 40, 120, and 190 MPa. The axial stress $\sigma_1$ is evaluated from the resulting force applied by the samples on the loading plates, and the volumetric strain $\epsilon_V$ is evaluated by monitoring the surface area of the domain delimited by the loading plates and the flexible membranes. A calibration by trials-and-errors provides rough optimum parameters for the contact law, which are provided in Table 1.

Figure 2 proposes a comparison between the experimental and numerical results. Predictions are satisfactory to the first order. The elastic response of Carrara marble is well represented, as well as its yielding beyond a certain stress level. The model also predicts, at least qualitatively, such key features as (a) the increase of the yielding stress with the confining stress, (b) the transition from a strength peak at low confining stresses to a work hardening at larger confining stresses, (c) a negative volumetric strain during the elastic part of the deformation followed by a dilatancy beyond the elastic domain, and (d) a strong reduction of this dilatancy as the confining stress
<caption>**Table 1**</caption>

Numerical and Mechanical Parameters of the Calibrated Discrete Model for Carrara Marble Simulation

<table><thead><tr><td>Parameter</td><td>Symbol</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>Contact stiffness</td><td>k</td><td>2.20 × 10<sup>14</sup></td><td>Pa/m</td></tr><tr><td>Ultimate strength</td><td>P<sub>lim</sub></td><td>50</td><td>MPa</td></tr><tr><td>Residual cohesion</td><td>P<sub>res</sub></td><td>10</td><td>MPa</td></tr><tr><td>Full damage gap</td><td>δ<sub>res</sub></td><td>4</td><td>μm</td></tr><tr><td>Contact loss gap</td><td>δ<sub>loss</sub></td><td>8</td><td>μm</td></tr><tr><td>Contact damping</td><td>γ</td><td>0.3</td><td></td></tr><tr><td>Marble density</td><td>ρ</td><td>2,700</td><td>kg/m<sup>3</sup></td></tr><tr><td>Confining stress</td><td>σ<sub>3</sub></td><td>5–40–120–190</td><td>MPa</td></tr><tr><td>Upper plate velocity</td><td>V<sub>y</sub></td><td>1</td><td>m/s</td></tr><tr><td>Axial strain rate</td><td>ε<sub>1</sub></td><td>26.2</td><td>1/s</td></tr><tr><td>Average grain size</td><td>d</td><td>180</td><td>μm</td></tr><tr><td>Number of grains</td><td></td><td>~17,000</td><td></td></tr><tr><td>Number of rigid degrees of freedom</td><td></td><td>~51,000</td><td></td></tr><tr><td>Number of deformable degrees of freedom</td><td></td><td>~2,000</td><td></td></tr><tr><td>Simulation time step</td><td>Δt</td><td>5</td><td>ns</td></tr></tbody></table>

increases. Some important behaviors are missing, though: under the present assumptions, the model is unable to predict the slow strength reduction beyond the stress peak at low $\sigma_3$ (and predicts a sharp stress drop instead), and the saturation of the yield stress beyond a certain level of confining stress. This last behavior is typical of the transition toward a pressure-independent plastic response of the marble under high confinement, while the proposed model sticks to a pressure-dependent Mohr-Coulomb-like model. Despite this discrepancy, however, the model is able to capture fairly well some features of the brittle-ductile transition in marble.

Figure 3 provides the evolution of damage distribution in the four samples during their loading history. It appears that the sample submitted to the lowest confining stress (5 MPa) fails in a sudden manner for an axial strain of ≈ 0.35%, following a pattern of two conjugate inclined and sharp fractures. The rest of the sample is virtually undamaged. At 40 Mpa, the failure is delayed to an axial strain of ≈ 0.50%, is less sudden, and occurs along two conjugate shear bands with a certain granular thickness, which concentrate almost all the damage of the sample. For the samples submitted to the two largest confining stresses, the failure patterns are different: they occur later (i.e., beyond an axial strain of ≈ 0.90%), and more progressively. They follow a more diffuse spatial distribution of the damage, which consists in a dense network of thin inclined shear cracks which eventually occupy the totality of the sample (except for the upper and lower areas which are stress-shadowed by the frictional contacts with the rigid loading plates). We thus recover with this model some key features observed in the experiments of Fredrich et al. (1989), namely a transition, as $\sigma_3$ increases, between a brittle and localized failure along fracture plane, and a more ductile and distributed failure.

Some close-in views of the microstructure are provided in Figure 4, as well as the corresponding damage distribution, for three situations: an intact sample, a sample that failed in a brittle manner, and a sample that failed in a ductile manner. The difference in the spatial distribution of the damage is evident, with localized open fractures separating intact regions for brittle failure, and dense arrays of inclined shear damage zones for ductile failure. Hence, despite the discrepancies observed in Figure 2 and commented above, the discrete model of marble presented in this section appears satisfactory enough for our purpose and is thus applied in the Laboratory Earthquake model of the next section.

### 3. Laboratory Earthquake Model

The numerical developments to be presented in the following sections are inspired by an experimental campaign reported in Aubry et al. (2020). These experiments are described in Figure 5a. The experimental samples are Carrara marble cylinders (40 mm diameter, 80 mm length), which are sawcut along a sliding interface oriented
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison between numerical and experimental (Fredrich et al., 1989) results in terms of differential stress ($\sigma_1 - \sigma_3$) and of volumetric strain $\epsilon_v$, at 5, 40, 120, and 190 MPa of confining pressure.</figcaption></figure>

by 30° with respect to the cylinders' axes (Figure 5b). These samples are placed within an autocompensated oil-loaded triaxial cell installed at the Laboratoire de Géologie de l'École Normale Supérieure, Paris, France. A confining stress $\sigma_3$ of 45, 90, or 180 MPa is applied to the sample, and a loading rate $V_y$ of 0.01 or 1 µm/s is applied vertically. Several quantities are monitored during the tests, including the axial stress $\sigma_1$ and the vertical shortening $\epsilon_1$ of the sample. Before experimental testing, the surfaces of the samples are prepared by flattening with a surface grinder and roughening either with #1200 or #50 wet silicon carbide grit sandpapers. This results in two surface configurations: “Smooth” surfaces (Figures 6a and 6b) with a roughness $S_q = 1.0$ µm, and “Rough” surfaces (Figures 6c and 6d) with a roughness $S_q = 12.7$ µm.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Failure patterns and damage distribution in the numerical reproduction of the experiments from Fredrich et al. (1989), for four different confining pressures (5, 40, 120, and 190 MPa) and at eight different axial strain levels—see also Movie S1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Detailed view of the numerical microstructure: (a) intact rock; (b) brittle and localized fracture of rock; (c) ductile and distributed damage of rock (damage color scale identical to Figure 3).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Reference experiment and model scale: (a) triaxial press used in Aubry et al. (2020); (b) geometry of the Carrara marble sawcut samples in Aubry et al. (2020); (c) 2D small fault model (S- fault, scale 1/10); (d) 2D medium fault model (M-fault, scale 1/4).</figcaption></figure>

The numerical samples inspired by this experiment are in 2D and are restricted to plane-strain kinematics. They represent the system at a reduced scale, namely at a 1/10 scale for the “Small fault” model (Figure 5c), and at a 1/4 scale for the “Medium fault” model (Figure 5d). They are presented in details in Figure 7. Each half-sample is composed of two areas with different mechanical treatments. They are essentially represented as an elastic continuum, discretized by a Multibody Meshfree approach which is very similar to the more classical Explicit Finite Element Method, albeit using Moving Least Square meshfree shape functions which improve accuracy and guaranty the continuity of the stress fields (Belytschko et al., 1994). An elastic constitutive law is applied to this continuum model, with the properties of Marble (Young's Modulus of 70 GPa and Poisson's ratio of 0.29, leading to a shear modulus of 27.1 GPa). A finite strain formulation is employed to handle the possibility of a non-negligible rotation of the elastic bodies during the vertical loading. In the first 200 µm of the subsurface of each half-sample, however, the numerical representation is switched to a discrete one, based on the model described and calibrated in the previous section. The only difference is a much smaller average grain size, close to 10 µm (Figure 7c), chosen as a very rough estimate of the average size of a large gouge particle (Passelègue, Schubnel, et al., 2016). Several parameters of the contact law (k, δ<sub>res</sub>, δ<sub>loss</sub>) are scaled with the grain size to keep an equivalent homogenized behavior as during calibration. A model with similar physics but using deformable grains showed a significant influence of the grain size on the material response (Falk et al., 2001). Numerical experiments on the calibration model (not reported here) however showed that the behaviors reported in Figure 2 are essentially independent from the grain size (although subjected to more variability if grains are larger), provided that this scaling is properly done. The discretization of the continuum zone is refined in the neighborhood of the discrete one, and a very highly cohesive law is applied between the two zones to ensure mechanical continuity (Figure 7b). The discrete surfaces to be put into contact are nominally flat, but are made rough by the presence of the irregular grains. As shown in Figure 6e, the resulting surface morphology is a compromise between the experimental rough and smooth surfaces, with a perfect macroscopic flatness (like the experimental smooth profile) and a large local roughness (like the experimental rough profile). The microstructure of Carrara marble is mostly based on calcite grains with a typical size of ≈ 150 µm (Figure 6f). It means that the discrete grains of the model are not intended to represent this microstructure, but rather to allow as fine a discretization as possible of the damaging, wear, and possible detachment of the matter at the rock surface during sliding events.

Two model sizes (Small and Medium) and three confining stresses (45, 90, and 180 MPa) are simulated, in order to isolate the effects of axial strain and of fault sliding distance. This results in six models called S45, S90, S180,
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Roughness and microstructure. (a) six illustrative profiles of 500 µm extracted from the smooth experimental samples from Aubry et al. (2020) (parallel and perpendicular to sliding direction); (b) topographic map of an initial smooth surface, with $S_q = 1.0$ µm (Aubry et al., 2020); (c) six illustrative profiles of 500 µm extracted from the rough experimental samples (parallel and perpendicular to sliding direction) from Aubry et al. (2020); (d) topographic map of an initial rough surface, with $S_q = 12.7$ µm (Aubry et al., 2020); (e) a typical zone of intact surface and microstructure of the numerical samples; (f) intact microstructure of marble before tests, with a typical grain size of 150 µm (Aubry et al., 2020). Scales for a, c, e, and f are identical.</figcaption></figure>

M45, M90, and M180. In the initial state, the two half-samples are not in contact. They are positioned between two elastic blocks (Figure 7a), which have two important roles: (a) introducing a compliance in the loading system, through an elastic modulus fitted to the experimentally measured apparatus stiffness, and (b) attenuating the propagation and reflexion of elastic waves through absorbing boundaries (Rajagopal et al., 2012). The elastic blocks are driven by rigid plates which control the vertical strain of the whole system. In the first computation stage, the contact status between the discrete grains is initialized to “intact,” and a confining stress $\sigma_3$ is applied on the lateral vertical walls of the half-samples and on the upper and lower rigid loading plates (with a correct- ing factor in order to yield an isotropic loading on the sample). The two half samples are therefore put into contact and stabilized in this isotropic state. The loading condition on the two rigid plates is then switched from pressure-controlled to displacement-controlled. The lower one is fixed in displacement, while the upper one is submitted to a constant vertical velocity $V_y$ chosen in order to maintain a vertical strain rate of $2.5 \ s^{-1}$, as measured between the rigid plates. When considering the steady-state shortening of the sample only, after the first elastic part, this strain rate is closer to $5.0 \ s^{-1}$ since the deformation in the loading system becomes almost stationary. This loading rate is orders of magnitude larger than the experimental one, but this artificial speed-up is rendered necessary by the high computational cost of the model. Loading rate is however low enough to be able to clearly distinguish quasi-static and dynamic episodes, as we elaborate later. The simulations are performed until an axial shortening of ~2.5% is reached in the samples at a final time of 5 ms, that is, after 25 millions of time steps of 0.2 ns. This corresponds to a final slip of ~200 µm for S-faults and ~500 µm for M-faults. Simulations were performed during several months on the computational cluster of the LaMCoS, Lyon, France. The numerical and mechanical parameters of the Laboratory Earthquake model are summarized in Table 2.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Main principles of the Laboratory Earthquake model (case of S-fault): (a) general view of the model, of its boundary conditions, and of its loading system; (b) zoom on the upper part of the model, including the elastic loading system in contact with the marble upper sample; (c) zoom on the fault zone, including the damageable surfaces after they were put into contact under isotropic confining stress. The triangular mesh is plotted to make the discretization fineness apparent, but does not correspond to the connectivity of the shape functions, which are meshfree.</figcaption></figure>

## 4. Effective Friction on the Fault

During the simulations, a large amount of data is collected at various rates, including the axial stress $\sigma_1$, the stress field and the velocity field in the sample, the position and the damage of the discrete grains, etc. To compute the effective friction $\mu$ (i.e., the instantaneous ratio between the average normal and tangential stresses on the fault), two different techniques are available. The first one is identical to the experimental one, and only relies on the knowledge of $\sigma_1$ and $\sigma_3$, and on the initial orientation angle $\theta = 60^\circ$ of the precut interface with the horizontal direction:

$$ \sigma_n = \left(\frac{\sigma_1 + \sigma_3}{2}\right) + \left(\frac{\sigma_1 - \sigma_3}{2}\right) \cos(2\theta) \quad (1) $$

$$ \tau = \left(\frac{\sigma_1 - \sigma_3}{2}\right) \sin(2\theta) \quad (2) $$

$$ \mu = \tau/\sigma_n \quad (3) $$

It is important to keep in mind that these formula postulate that (a) the angle $\theta$ does not change during loading, and that (b) the axial and confining stresses are the principal components of the stress tensor, at least when it is averaged on the whole system. Another approach consists in evaluating the normal and tangential stresses on each point of the fault surface. This is obviously out of reach in experiments, but can be well approached in the present model by considering the planar frontier between the continuous and discrete parts, for example, in the lower half-sample. At any node of this frontier, the complete stress field is known through the equations of plane strain elasticity implemented in the continuous part, and one can therefore write (following the frame notations of Figure 7, where $\sigma_{nn}$ is the normal stress on the fault plane and $\sigma_{nt}$ is the tangential stress):

$$ \sigma_{nn} = \sigma_{xx} \cdot \sin^2 \theta + \sigma_{yy} \cdot \cos^2 \theta - 2\sigma_{xy} \cdot \cos \theta \cdot \sin \theta \quad (4) $$
<table><caption>**Table 2**

Numerical and Mechanical Parameters of the Laboratory Earthquake Models</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>S-faults</th><th>M-faults</th><th>Unit</th></tr></thead><tbody><tr><td>Contact stiffness</td><td>k</td><td>3.92 × 10<sup>15</sup></td><td></td><td>Pa/m</td></tr><tr><td>Ultimate stress</td><td>P<sub>lim</sub></td><td>50</td><td></td><td>MPa</td></tr><tr><td>Residual cohesion</td><td>P<sub>res</sub></td><td>10</td><td></td><td>Mpa</td></tr><tr><td>Full damage gap</td><td>δ<sub>res</sub></td><td>0.2</td><td></td><td>μm</td></tr><tr><td>Contact loss gap</td><td>δ<sub>loss</sub></td><td>0.4</td><td></td><td>μm</td></tr><tr><td>Contact damping</td><td>γ</td><td>0.3</td><td></td><td></td></tr><tr><td>Marble density</td><td>ρ</td><td>2,700</td><td></td><td>kg/m<sup>3</sup></td></tr><tr><td>Scale</td><td></td><td>1/10</td><td>1/4</td><td></td></tr><tr><td>Height</td><td>H</td><td>8</td><td>20</td><td>Mm</td></tr><tr><td>Width</td><td>W</td><td>4</td><td>10</td><td>Mm</td></tr><tr><td>Confining stress</td><td>σ<sub>3</sub></td><td>45–90–180</td><td></td><td>Mpa</td></tr><tr><td>Upper plate velocity</td><td>V<sub>y</sub></td><td>40</td><td>100</td><td>mm/s</td></tr><tr><td>Axial strain rate</td><td>ε<sub>1</sub></td><td>~5.0</td><td></td><td>1/s</td></tr><tr><td>Average grain size</td><td>d</td><td>10</td><td></td><td>μm</td></tr><tr><td>Fault length</td><td></td><td>~6</td><td>~16</td><td>mm</td></tr><tr><td>Number of grains</td><td></td><td>~31,000</td><td>~83,000</td><td></td></tr><tr><td>Number of rigid Degrees of Freedom</td><td></td><td>~93,000</td><td>~249,000</td><td></td></tr><tr><td>Number of deformable Degrees of Freedom</td><td></td><td>~37,000</td><td>~115,000</td><td></td></tr><tr><td>Simulation time step</td><td>Δt</td><td>0.2</td><td></td><td>ns</td></tr><tr><td>Number of time steps</td><td></td><td>2.5 × 10<sup>7</sup></td><td></td><td></td></tr></tbody></table>

$$ \sigma_{nt} = (\sigma_{xx} - \sigma_{yy}) \cdot \cos \theta \cdot \sin \theta + \sigma_{xy} (\sin^2 \theta - \cos^2 \theta) \quad (5) $$

$$ \mu = \int \sigma_{nn} / \int \sigma_{nt} \quad (6) $$

In these expressions, the orientation $\theta$ of the fault is monitored and updated during simulation. The integrals in Equation 6 are on the whole fault surface. The Equations 3 and 6 are used to plot two evaluations of the effective friction for the S180 simulation in Figure 8. We observe typical patterns of stick-slip, which are characteristic of the seismic cycle in the case of highly coupled faults, both in the field and in the lab. Stress drops are of various sizes, and seem consistently described by both techniques. However, while the orders of magnitude and local variations are consistent between the two approaches, there is a general drift of the first method with respect to the second. This drift is to be attributed to finite strain effects (i.e., the change of orientation of the fault as the sample is compressed) and to the development of a non-zero $\sigma_{xy}$ term. This shear stress develops because of the use of a frictional interface between the sample and the loading system (in agreement with the experimental setup), which restrains the relative horizontal motion of the two half samples as slip develops on the fault (see Figure 11d, commented in Section 5). As compression develops, the axial and confining stresses $\sigma_1$ and $\sigma_3$ are no longer the principal components of the average stress field in the sample. Hence, while the simpler experimental approach can readily be used to analyze stress drop patterns, it is slightly misleading when it comes to the average fault friction.

The mechanical results of the whole numerical campaign are presented in Figure 9, where effective friction (evaluated based on the average normal and tangential stresses on the fault, Equation 6) is plotted against axial strain (evaluated from the motion of the upper and lower boundaries of the sample) for the cases S45, S90, S180, M45, M90, and M180. In all cases, samples first experience an elastic loading until a first friction peak (between effective frictions of 0.3 and 0.4) and a large stress drop, followed by a stick-slip regime around values of 0.2–0.3. A very large number of events of various sizes can be observed in the insert of Figure 9 for both model sizes and for the three confining stresses.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Effective friction on the S180 fault as monitored (i) through an estimation based on the axial stress $σ_1$ applied by the loading system, and (ii) through an integration of the stresses measured on the continuous border of the lower half-sample.</figcaption></figure>

These results are summarized in Figure 10 and compared with the experimental results from Aubry et al. (2020). Only the experimental samples compressed at the larger loading rate (1 µm/s) are considered, since slower loadings reported in Aubry et al. (2020) lead to different behaviors, most likely due to healing effects during loading. Since no healing is implemented in the model and since the numerical loading rate is orders of magnitude above the experimental ones, we thus discarded the experiments at the lower loading rate of 0.01 µm/s for the sake of fair comparison. Figures 10a and 10b provide the experimental curves of effective friction against axial strain, for smooth and rough faults respectively (see illustrative roughness profiles in Figure 6). As reported in Aubry et al. (2020), the results for smooth and rough faults are very different. On the smooth faults, the tests

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Effective friction evolution with axial strain in the six numerical samples.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparison between experiments and simulations. (a) experimental effective friction for smooth faults from Aubry et al. (2020) (tests were doubled for $\sigma_3 = 45$ and 90 MPa); (b) experimental effective friction for rough faults from Aubry et al. (2020); (c) comparison between the experimental and numerical average post-peak friction (experimental values averaged on the bold parts of curves in a and b); (d) comparison between the experimental and numerical effective friction at the first peak (note that no peak were observed in the experiments with rough faults).</figcaption></figure>

with confining stresses of 45 and 90 MPa lead to a peak effective friction followed by a large stress drop, and by either a stable sliding or a succession of small and regular stick-slip events. At 180 MPa, there is also a first large stress drop, but it is followed by a second one and by a continuous increase of the effective friction up to ~0.5. In the case of rough faults, for the three confining stresses, there is no friction peak and effective friction increases steadily until the end of the tests, reaching ~0.8 in the case $\sigma_3 = 45$ MPa. For all these simulations, an average effective friction was approximately computed by considering “steady state” sections of the curves (chosen rather arbitrarily in the areas where friction tends to stabilize, and plotted in bold lines in Figures 10a and 10b). Peak effective frictions were also extracted for the smooth fault cases. These results are compared with the numerical ones in Figures 10c and 10d respectively. Figure 10c shows that the numerical average friction slightly decreases with the confining stress. It is also slightly larger for the S-faults than for the M-faults. The numerical values are much lower than the experimental ones obtained on the rough fault, but are consistent with the experimental results if we consider the smooth fault for $\sigma_3 = 45$ and 90 MPa. There is a large discrepancy with the experiment at 180 MPa, though, which indicates that the numerical model is unable to represent the change of behavior from brittle to ductile at large confining stresses and the onset of bulk plasticity. This was to be expected because of the elastic constitutive law implemented in the continuous regions of the model. Peak effective frictions (Figure 10d) follow the same trend, with numerical values in good agreement with experimental ones except for the case $\sigma_3 = 180$ MPa, for the same reason as before. We also note that the first peak values are systematically larger in the numerical S-faults than in the M-faults.

## 5. Heterogeneity and Sliding Events

As shown in Figures 8 and 9, the numerical models accommodate sliding by a succession of slip events of various sizes. Figure 11 provides some selected views of the stress field in the sample at several stages of the simulation M180. The stress state is represented through the ratio between the tangential stress $\sigma_{nt}$ and the normal stress $\sigma_{nn}$ oriented along the direction of the fault plane, as provided by Equations 5 and 4 respectively. On the fault plane itself (i.e., at the frontier between the continuous and discrete regions of a given half-sample), this ratio is equal to the effective friction as provided by Equation 6. Figure 11a shows that, in the first stage of the simulation, before the beginning of the axial loading, the stress ratio is close to zero in most of the sample. This is in good agreement with the initially isotropic stress state created by the confining stress $\sigma_3$ alone. At the end of the elastic loading (Figure 11b), the stress ratio has increased to an average value of ~0.3. We observe however that it is highly heterogeneous in the neighborhood of the fault plane. Since no large slip event has yet occurred at that stage, and in the absence of roughness, this initial stress heterogeneity cannot be attributed to the sliding history or to load concentration around asperities. It must therefore be related to an intrinsic heterogeneity in the surface response
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Evolution of the stress field in the sample M180 at several loading stages. Stress state is plotted as the ratio between the tangential stress $\sigma_{nt}$ (Equation 5) and the normal stress $\sigma_{nn}$ (Equation 4) oriented along the fault frame {$n, t$} (Figure 7). In the neighborhood of the fault plane, this is equivalent to an effective friction. (a) before axial loading; (b) just before the first stress peak; (c) just after the first stress peak; (d) during the steady-state stick-slip regime; (e) evolution of average effective stress in time, and time marks for figures (a) to (d), see also Movie S2.</figcaption></figure>

to tangential loading, probably because of the random character of the polygonal microstructure introduced in the degradable part of the half-samples. After the first major sliding event (Figure 11c), this heterogeneity is still present, with a similar level of variability. The drop in stress ratio is clearly visible during this sliding event. As slip develops on the interface and a large number of small and large events occur (Figure 11d), we observe two interesting evolutions. First, the stress ratio profile on the fault seemingly follows slightly different patterns. Heterogeneities seem both more intense, less evenly distributed, and with a larger characteristic spatial period, with fault areas of several hundreds of µm with very low (>0.1) or very high (>0.5) stress ratios. It means that the fault strength, while averaging at a stress ratio close to 0.25, is extremely variable along the fault. Some patches may play the role of asperities (i.e., locus of larger shear resistance with respect to the fault average) without the need of an initial roughness. Second, we observe a breakdown of the large-scale homogeneity and symmetry of the stress field in the sample, with high shear stresses in the upper-left and lower-right corners. This observation is related to finite-strain effects, because the horizontal motion related to the slip between the two-half samples is restrained by the frictional contact at their interface with the elastic loading blocks above and below. This observation explains the discrepancy observed in Figure 8, and likely occurs in the lab as well since the interfaces between the sample and the loading system are not lubricated.

The heterogeneities observed in the stress field in Figure 11d are the product of (a) granular phenomena within the interface (illustrated in the next section) and of (b) a complex sliding history. In Figure 12, we analyze three sliding events in the simulation M180. For each event, we provide (from left to right) a view of the differential
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Three examples of sliding events of varying sizes in the case M180. (a) large event on the whole fault with, from left to right (i) differential Coulomb stress between before and after the event, (ii) zoomed view on a selected area of the fault, and (iii) axial stress measured in the loading system with black dots indicating the times where stress fields were extracted for differential Coulomb stress calculation; (b) intermediate-size event in the lower half of the fault; (c) small event completely localized inside the fault; (d) axial stress history of the M180 case with marks for events (a) to (c), see also Movie S3.</figcaption></figure>

Coulomb Stress in the sample before and after the event, a close-in view of a selected part of the fault, and the history of the axial stress $σ_1$ monitored in the loading system in the temporal vicinity of the event. The instants at which the stress field is sampled to evaluate the differential Coulomb stress are marked by black dots. The first event (Figure 12a) corresponds to the main stress drop at the end of the elastic loading, for which the entire fault slipped for about 25 µm. This event appears very clearly on the $σ_1$ plot, with an axial stress drop close to 36 MPa. We observe a general and rather uniform decrease of the Coulomb stress in the whole sample (by about 15 MPa), but a large heterogeneity on the fault plane, with local variations ranging from -50 to +30 MPa on fault
patches with a length of ~200–400 µm length. Despite the general stress drop and the large sliding distance, it is striking to observe that some local areas of the fault are actually much more stressed after the event than before. We therefore have a first source of stress heterogeneity on the fault, which spontaneously occurs during single large sliding events.

Figure 12b shows an intermediate-size sliding event, which roughly ruptured one third of the fault and reached its lower boundary. This rupture is visible on the axial stress plot, with a stress drop of about 1 MPa as measured in the loading system. Between the two sampling times used to compute the field of differential Coulomb stress, however, the value of $σ_1$ increased because of the continuous loading. This is the reason why the differential Coulomb stress field exhibits a general stress increase in the sample, close to ~2 MPa. In the area of the fault which unlocked, however, the Coulomb stress drop is very clear, and locally reaches more than 10 MPa. Just above the ruptured area, we observe a clear evidence of stress concentration on the fault, with a Coulomb stress increase of about 10 MPa. Two different scales of heterogeneity are visible, since the patterns observed in Figure 12a are still present and are superimposed with the lower-wavelength heterogeneity arising from the partial failure of the fault. Finally, the event of Figure 12c is very small. It cannot be detected based solely on the $σ_1$ history, since it is too local and only concerns about 10%–15% of the fault length. In that case, we observe that the Coulomb stress drop in the ruptured zone is rather homogeneous, with an average drop of 2 MPa. There is still some evidence of stress concentration at the two extremities of the sliding region. Around the ruptured patch, we observe that the fault remained locked and the Coulomb stress kept increasing during the sampling interval. From this observation we may conclude that the stress heterogeneities observed in Figure 11 are both related to the presence of partial sliding on the fault and to spontaneous spatial variability of the stress drop during large events.

Figure 13 proposes a more detailed view of the kinematics of the event described in Figures 11b, 11c, and 12a. The instantaneous velocity field in the sample is first corrected by subtracting the affine velocity field related to the continuous axial loading of the sample, and then projected in the direction tangential to the fault plane. This field is plotted at five selected times in Figures 13a–13e. In Figures 13a, 50 µs before the mainshock, we observe moderate levels of velocity in the sample, related to elastic waves of previous minor events which have not completely dissipated. Despite this non-null velocity, it is clear that the interface is completely locked since there is a perfect continuity of the velocity field across it. In Figure 13b, about 3 µs before the mainshock, this perfect coupling no longer holds. We observe localized slips in the lower and upper areas of the fault, with some patches which remain locked in between. Since the stress peak is still to come, these localized slips can certainly be assimilated to foreshocks. In Figure 13c, we observe the mainshock, with a sliding velocity (i.e., a velocity jump between the two half samples) of the order of 10 m/s and a complete sliding of the interface. As shown in the zoomed view, the velocity jump occurs within the granular part of the fault (see next section). It is extremely localized on some areas (upper part of the insert of Figure 13c) and more diffuse on some others (lower part of the same insert). About 4 µs after the end of the mainshock (Figure 13d), the fault is mostly locked again except in some localized areas where small and short slip events still occur (see insert in Figure 13d), which we can probably qualify as aftershocks. Finally, 40 µs after the mainshock (Figure 13e), the fault is completely locked, there is no discontinuity of the velocity field across it, and the remaining non-null velocity is related to the elastic waves produced by the aftershocks still bouncing in the sample and its loading system before complete attenuation. This figure therefore draws the picture of very complex phenomena which are reminiscent of typical features of the seismic cycle as observed in the field: a locked fault, a typical sequence foreshocks-mainshock-aftershocks with large heterogeneities in space and time, and a relocking of the fault. As shown in Figure 9, this cycle occurs a very large number of times in each simulation. It seems, however, that many events deviate from this general scheme and occur without foreshocks or aftershocks, but this will have to be confirmed by more detailed studies in the future.

## 6. Damage and Shear Accommodation

Since the elastic regions of the model are perfectly homogeneous, a large part of the spatial and temporal complexity observed in the previous section is to be attributed to the discrete regions of the sample. These regions are statistically homogeneous (meaning that the average size, connectivity, orientation, etc. of the polygonal grains is constant on the whole fault) and are nominally flat (apart from the grain-scale roughness, Figure 6), but there is nevertheless a randomness in the geometrical arrangement of the grains, which is sufficient to trigger the apparently chaotic behavior described in the previous section. In this section, we therefore focus on the granular
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Velocity fields for the large sliding event of Figure 12a in simulation M180. Velocities are corrected from the affine component related to the axial loading, projected in the tangential direction to the fault, and plotted with a logarithmic color scale. Subplots (a–e) correspond to different times, marked on the axial stress history in (f)—see also Movie S3.</figcaption></figure>

phenomena which occur within the interface, and which primarily control the rheology and therefore the frictional response of the fault. An important phenomenon to observe is the progressive damaging of the interface and the progressive build-up of a layer of granular gouge. Figure 14 shows the evolution in time of the damage at selected locations of the M45, M90, and M180 faults. In the M45 test, the distribution of damage in the first stage of the test (i.e., for $t = 0.25–0.50$ ms) is heterogeneous, with damaged patches extending a few tens of µm, separated by small intact region. As sliding develops, however, we observe that the level of damage tends to become more homogeneous along the fault, with the development of a fully damaged granular layer in the interface. At $t = 5$ ms, this layer presents a homogeneous thickness of ~80 µm, and is bounded by almost intact discrete regions on both sides. The simulation M180 shows slightly different patterns, with the appearance of a number of fractures which originate at the fault and propagate up to the end of the discrete region in a sub-perpendicular direction. These fractures come in addition to the development of the damaged granular layer, which starts later than in the M45 case (because the first major slip event only occurs around $t = 0.96$ ms in the present case) but develops faster. At $t = 5$ ms, it has an average thickness of ~120 µm. In contrast with the M45 case, this thickness is not homogeneous along the fault, and the central layer is bounded by highly damaged regions on both sides,
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Evolution in time of the damage patterns on selected discrete regions of the faults M45, M90, and M180. Numbers in white indicate time in ms—see also Movie S4.</figcaption></figure>

with patterns similar to those of Figure 4c (obtained under a similar confining stress). The M90 case seems to adopt an intermediate behavior between M45 and M180.

In Figure 15, the profile of relative damage across the interface is averaged along the whole fault for each stored state of the model S180, and these profiles are then plotted along time (the relative damage being color-coded from blue to red). On the same figure, the curve of the average tangential stress on the fault (Equation 2) is plotted on the same time axis for the sake of comparison. We observe the progressive development of the damage in the interface during the simulation, and the growth of a central layer with a complete damage, for which all cohesive bonds between grains were broken. It is evident from this plot that most of the damage takes place during the first large sliding events that occur on the fault (e.g., at times close to 1.10, 1.20, and 1.35 ms). Two important points nevertheless deserve to be noticed. First, there is a clear and steady increase of the relative damage during the

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Evolution of the average damage profile of the interface of the simulation S180, superimposed with the history of the average tangential stress on the fault.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Evolution of the average damage profiles of the simulations S45, S90, S180, M45, M90, and M180.</figcaption></figure>

initial phase of elastic loading, where only minor events occur. And second, after a time of ~2 ms, large events keep occurring on the fault but do not result any more in sudden increases of the relative damage. As illustrated in Figure 15, we can therefore distinguish three different stages regarding fault damage. (a) The elastic stage, where there is a steady increase of damage at a moderate rate, sometimes accelerated by minor foreshocks, (b) a transient stage where several large events lead to sudden increases of the damage and to the appearance of a central gouge layer, and (c) a mature stage where seismic cycles keep occurring but damage develops at a much slower and more regular rate. Similar evolutions are observed for the other models (Figure 16).

The evolution of damage is quantitatively analyzed in Figure 17. Figure 17a presents five average damage profiles for the M90 fault, at different times of the simulation. At $t = 0.5$ ms, the first large sliding event has not occurred yet, and the damage profile across the fault is Gaussian. However, from $t = 1$ ms to the end of the simulation, we observe a central region of full damage surrounded by two Gaussian-like damage decays. The central fully damaged layer thickens progressively as the simulation progresses. We define a characteristic Damage Thickness $T_D$ as the interface thickness for which the relative damage is larger than 0.4 (Figure 17a), and follow the evolution of this parameter in time in Figure 17c for the six simulations. The six curves follow the same evolution, with an initial stage of very limited damage, a sudden and large increase of the damage thickness, and a progressive shift toward a linear increase of this stiffness with simulation time. This is in good accordance with the observations of Figure 15. Figure 17c also shows that damage thickness is larger for simulations with large confining stresses. For a given value of $\sigma_3$, it also tends to be slightly larger for the M-faults than for the S-faults, although the evolutions are similar for the two model sizes. This last observation is particularly interesting, because at a given time (i.e., for a given value of the axial strain $\epsilon_1$), the sliding distance is about 2.5 times larger for M-faults than for S-faults. It therefore appears that the damage thickness is not related to slip, but to strain. To validate this observation, we extracted the slopes of the linear regions of the six curves
<figure><img src="image_19.png" /><figcaption>**Figure 17.** A quantitative view of damage and shear accommodation. (a) average damage profiles at different times of the simulation M90; (b) average shear rate profiles at different times of the simulation M90; (c) evolution of the Damage thickness $T_D$ (corresponding to a damage of 0.4, see a) in time for the six simulations; (d) evolution of the Shear Thickness $T_S$ (corresponding to 20% of the maximum shear rate, see b) in time for the six simulations; (e) rate of increase of $T_D$ with the axial strain, as measured by the linear fits of (c); (f) rate of increase of $T_S$ with the sliding distance, as measured by the linear fits of (d).</figcaption></figure>

of Figure 17c (corresponding to the “mature stage” of Figure 15), and used them to compute a damage rate $\dot{T}_D$, expressed as:

$$ \dot{T}_D = \frac{\delta T_D}{\delta \varepsilon_1} \qquad (7) $$

Where the evolution of $\varepsilon_1$ is averaged on the same time period. $\dot{T}_D$ is expressed in µm/µdef, where a µdef is a strain unit (µm/m). The results are plotted in Figure 17e, and confirm that, in the stabilized regime of the faults, the damage thickness rate is proportional to the confining stress $\sigma_3$. An increment of damage thickness $\delta T_D$ can thus be written as:

$$ \delta T_D \propto \delta \varepsilon_1 \cdot \sigma_3 = \delta W \qquad (8) $$

where $\delta W$ is an increment of elastic energy density stored in the sample during loading. Damage evolution is therefore essentially controlled by a constant share of the energy budget of the system, regardless of the size of the fault and of the confining stress.

Figure 17b shows profiles of shear rate in the fault interface at several times of the simulation M90. Since the instantaneous shear rate is close to zero when the fault is locked and has very large values when it unlocks, and since it varies in space because of the complexity of the sliding history (Figure 12), the profiles plotted in Figure 17b result from a time-averaging of the grains velocities on a sufficient time window. We observe a Gaussian-like distribution of the average shear rate, with a central maximum value which tends to decrease as sliding develops, while the profile gets wider. We define a characteristic shear thickness $T_S$, for which the average shear rate is larger than 0.2 times its maximum value at the center of the interface. The quantity $T_S$ is plotted over
time in Figure 17d for each of the six simulations. In contrast with the damage thickness $T_D$, the evolution of $T_S$ is linear in time during the whole simulation (early instants are discarded because sliding has not occurred yet). Shear thickness appears to be largely independent on the confining stress, but to be larger and to increase faster in time for the M-faults as compared to the S-faults. This size-dependence is certainly related to the fact that, in average, the M-faults slide about 2.5 times faster than the S-faults (for an equivalent strain rate as applied by the loading system). To validate this hypothesis, we fit a linear equation on the evolution of $T_S$ and we define a rate of increase of the shear thickness per sliding distance $\dot{T}_S$:

$$ \dot{T}_S = \frac{\delta T_S}{\delta U} \qquad (9) $$

where $U$ is the relative tangential motion between the two half samples (averaged on the same time interval as $T_S$). The results are summarized in Figure 17f, which unambiguously shows that the rate of increase of the shear thickness is both independent from the sample size and from the confining stress, and depends solely on the sliding distance. While the condition $T_D > T_S$ should necessary be fulfilled, the model therefore indicates that shearing is not accommodated on the whole damage thickness, and that a large part of the damaged area remains linked to the host rock and does not participate to sliding. This is particularly the case when the confining stress is large compared to the rock strength.

# 7. Discussion

## 7.1. A Numerical View of the Seismic Cycle

The numerical results provided by the laboratory earthquake model bring a clear view of a succession of events that can probably be considered as a good representation of the seismic cycle in the lab. As axial loading starts, a tangential stress $\sigma_{nt}$ builds up in a rather homogeneous way in the sample (Figure 11), with some heterogeneities in the vicinity of the fault plane. This tangential stress is endured by the fault interface, although some progressive damaging starts to occur well-before the first major sliding event (Figures 14 and 15). No major motion occurs until $\sigma_{nt}$ reaches the strength of the fault on some localized patches, which yield and trigger minor slip events (Figure 13b). Shortly after these foreshocks, a first mainshock occurs on the whole interface. In the M180 sample, it is characterized by a maximum sliding velocity of ~10 m/s, a sliding distance of ~25 µm, and a duration of ~10 µs. This mainshock corresponds to a drop in tangential stress of about 30 MPa as averaged on the whole fault interface, although this drop is highly spatially variable (Figure 11c). This results in a new stress state in the sample, with a general stress level lower than before the event but with a similar heterogeneity (Figure 12a). The sliding event is concomitant with a sudden and large increase of the damaging of the interface, and with the creation of a central layer of fully damaged rock (fault gouge, Figure 14). The mainshock is followed by a small number of minor events, which might qualify as aftershocks (Figure 13d), until the fault locks again (Figure 13e). The loading cycle then resumes, with a succession of a large number of sliding events of various sizes, some of which concerning the whole fault, but most of which concerning only a subset of it (Figure 12). During the first few large events, we observe sudden increases of the damage thickness within the interface (Figure 15), but the rate of thickening of the damage zone then decreases and seems to stabilize to a constant value (Figure 17c). Meanwhile, the accumulation of wear material in the interface leads to the creation of a tribological third body, which accommodates the shearing between the two pieces of solid rock (Figure 17b). We observe that this accommodation thickness increases as the seismic cycles develop (Figure 17d), but does not follow the same rate as the interface damage. Under the assumptions of the present model, no subsequent event is as large as the first one. After a large number of cycles and a total sliding distance close to 500 µm (for the M-faults), the stress state in the sample is still very heterogeneous close to the fault, and the intensity and the spatial wavelength of these heterogeneity patterns seem to increase with the total slip distance (Figure 11d).

## 7.2. Comparison With the Reference Experiment

One interesting point is that we did not postulate any ad hoc weakening friction law in order for the system to respond to loading through stick-slip patterns. All the mechanical parameters of the model arise from an independent calibration on bulk marble, and the interface is then left to evolve spontaneously toward its preferred sliding behavior. In the present case, stick-slip patterns arise naturally from the interactions between rock elasticity
and (a) its calibrated fracture properties (i.e., asperity breakage in the very first stages of the loading), (b) its progressive damage, (c) the production of a gouge in the interface, and (d) the granular rheology of this third body layer (Aharonov & Sparks, 2004; Houdoux et al., 2021; Lherminier et al., 2019), all of which being consequences of the bulk properties of the host rock. This observation confirms that the typical friction laws used in numerical simulations of earthquake nucleation and propagation (e.g., slip-weakening or rate-and-state friction models) are mere consequences of a large number of local phenomena taking place within the interface (Casas et al., 2022; Mollon et al., 2021), and are by no means to be considered as first-principles.

It is worth mentioning, that the model fails at that point to reproduce accurately some of the experimental results of Aubry et al. (2020). Some comparisons with the experiments of reference are satisfactory: this includes the values of the peak and average friction (at least when the comparison is made with the case of “smooth” surfaces and for confining stresses of 45 or 90 MPa), which are nicely reproduced by the model (Figure 10). Stick slip patterns are also present in the experiments in those cases, albeit with much smaller stress drop values. The model however fails to reproduce experiments in the case of “rough” surfaces, and in the case of a confining stress of 180 MPa. Besides, the fact that large confining stresses lead to larger damage thickness is evident under the assumptions of our model, and is in line with the general spirit of the Archard law, but is much less evident in experiments because large confining stress switch the major accommodation site from the interface to the bulk marble (Aubry et al., 2020). The major reasons for these discrepancies are easy to identify: (a) the surface roughness of the model is neither identical to that of the “rough” nor of the “smooth” experimental surfaces, and is somewhere in-between (Figure 6), (b) there is no healing or aging of the interface during the simulation, and (c) the constitutive law of the continuous part of the model is purely elastic while bulk plasticity of marble is expected to take over at large confining stresses (Figure 2). The necessary improvements shall be introduced in future versions of this model.

One important question is that of the grain size. In our model, the gouge grains are inherited from the discrete representation of the intact marble surfaces (Figure 6e). They are very homogeneous in size, with an equivalent diameter close to 10 µm. This happens to be close to a typical gouge grain size reported in Passelègue, Schubnel, et al. (2016), but cannot claim to reproduce in a realistic manner the typical self-affine granulometry of fault gouges, which generally include large populations of much smaller grains. Since sliding dynamics predicted by the model result from grain-scale processes, we may wonder to what extent this assumption affects the results. Another type of grains that were not properly reproduced in the model are the mineral grains composing the microstructure of the intact rock (Figure 5f). A future model could include such features, in relation with the concepts of thermal cracking analyzed in Passelègue, Spagnuolo, et al. (2016). As demonstrated experimentally in Aubry et al. (2018), heat creation in the interface during the coseismic sliding has important consequences on its frictional response. This may include pervasive melting of the gouge layer (Mollon et al., 2021) and subsequent quenching and welding of the interface (Aubry et al., 2018), but also off-fault damage related to heat diffusion: thermal expansion anisotropy and heterogeneities among mineral grains of the bulk rock may lead to large stress concentrations when a large amount of heat diffuses away from the interface, and the resulting damage and surface creation are far from negligible (Passelègue, Spagnuolo, et al., 2016). This feature could be tackled using special refinements of DEM including heat diffusion and thermal expansion of grains.

### 7.3. Wear and Damage

As demonstrated experimentally in Aubry et al. (2018), fault structure evolves along the seismic cycles. One reason for this evolution is the occurrence of wear. Wear rate is a notoriously complicated quantity to measure, either in the lab or in the field. One of the reasons is probably that this quantity is rather ill-defined, as illustrated in the present study. Wear rate could be defined based (a) on the amount of matter that is completely damaged (but this concept is difficult to apply rigorously outside of an idealized rock damage model in numerical modeling), (b) on the quantity of gouge material participating to the shear accommodation in the interface (but this requires kinematic measurements), or even, in a tribological framework, as the quantity of matter that is definitely ejected from the contact (Godet, 1984) (and it would therefore be equal to zero in our model and in real faults). It is generally recognized in tribological practice that wear occurs in two phases: running-in and steady-state. Boneh et al. (2014) focused their work on that aspect, using a rotary shear apparatus on various rock types with normal stresses up to 7 MPa (but generally close to 1–2 MPa) and sliding velocities up to ~1 m/s (but generally close to 0.01–0.1 m/s). Among other quantities, they evaluated wear by measuring fault-normal displacement
during sliding. They noticed a first intense wear stage corresponding to asperity breakage and flattening, in accordance with the classical view of wear proposed by Archard (1953), postulating that the initial contact area is only a small fraction of the nominal one (Bowden & Tabor, 1950). These processes do not occur in our model because the initial surfaces are nominally flat and the normal stress is way larger. In Boneh et al. (2014), this very short stage was followed by a running-in period, that is, a transition period of 1–3 m of sliding, during which both friction and wear rate decreased until a final stage of steady-state, where the wear rate was found roughly constant. These wear kinetics, already proposed experimentally in W. Wang & Scholz (1994) and numerically in Guo & Morgan (2007), are in good agreement with our numerical predictions (Figure 17). The running-in period was interpreted in Boneh et al. (2014) as the transition from a two-body to a three-body sliding regime through the production of an established layer of fault gouge. In our simulation, it would correspond to the fault transitioning from an initial state where the only way to initiate sliding is to break the geometrical interlocking of the grain-scale roughness, to a mature state where sliding can occur by shearing of a granular layer, as interpreted in Aubry et al. (2018). This explains why the first peak (essentially controlled by the ultimate strength $P_{lim}$ of the intergranular bonds) is systematically larger than the following ones (essentially controlled by the shear strength of the gouge layer). As demonstrated numerically in Fillot et al. (2007), the stable third-body thickness in an open tribosystem corresponds to a balance between a source flow (i.e., a flow of matter detached from the surfaces and released in the interface) and an ejection flow (i.e., a flow of matter definitely lost to the system). This thickness (and therefore the wear rate) is thus both a property of the materials and of the system, because it depends on the time needed for third body to reach the contact boundaries and to get ejected. These concepts readily apply to the experiments of Boneh et al. (2014). In a closed system (as in the present model and in real faults), however, such ejection is impossible and the gouge thickness keeps increasing. Tribological experiments show that the source flow is expected to decrease in time as the larger third body thickness protects the surfaces by enduring most of the shearing (Hsia et al., 2020). This process is similar to shear localization in dense granular materials, where shear naturally occurs in a layer with a certain material-dependant thickness, leaving the rest of the system relatively undisturbed (Mollon et al., 2021). This is also what we observe in the present study. Hence, in our case, the running-in period is not a transition between a two-body and a three-body tribosystem (since a small amount of third body is pervasively present in the interface from the first major slip event, Figures 14 and 15), but it corresponds to the time needed to build a sufficiently thick layer to protect the surfaces. This layer is evaluated to 10–20 grains in our discrete model (about 100–200 µm in Figure 17c), but extrapolation of this estimate to lab or field situations would certainly be premature.

Our quantitative estimates of the steady-state wear rate (as extracted from Figure 17f) is close to ~40,000 µm/m (meaning that for 1 m of sliding, active gouge thickness would increase by 40,000 µm). This estimate is based on the sheared thickness of the gouge layer, which is only one imperfect way to quantify wear. During their experimental campaign, (Boneh et al., 2014) reported very various values of steady-state wear rates, ranging from 0.4 to 961 µm/m, with an average value close to 60 µm/m. This is about three orders of magnitude smaller than our predictions, but their tests were run under normal stresses close to 1–2 MPa. Figure 17f tends to indicate that the $T_s$ value is independent from the confining stress, but this contradicts Figure 17e which shows that the damaged thickness vanishes with the confining stress. Hence, the extrapolation of Figure 17f toward confining pressures of the order of those reported experimentally in Boneh et al. (2014) is problematic. We can only guess that, if no damage occurs, no gouge thickness will develop, meaning that $T_s$ should vanish at low confining stress just as $T_D$ does. Applying a simple Archard law (i.e., wear rate would be roughly proportional to normal stress on the contact), we find that the experimental and numerical estimates are not completely inconsistent. We can also speculate that our numerical wear rate is not independent from the chosen grain size in the degradable layer. While this grain size (~10 µm) is not unrealistic, it remains rather arbitrary, and its implication on the numerical wear rate will have to be determined in future numerical experiments.

In the experiments used as a reference for our model (Aubry et al., 2020), the wear rate was not quantified but *post-mortem* observations allowed to determine that the same amount of variability as in Boneh et al. (2014) should be expected, depending mostly on the initial surface roughness: smooth faults led to a final aspect that was either mirror-like with striations (for confining stresses of 45 and 90 MPa) or matte without striation (180 MPa), without much gouge produced; rough faults led to mirror-like surfaces with a thick layer of gouge under moderate confining stresses (45 and 90 MPa), and to matte surfaces without gouge in the case of 180 MPa. The absence of gouge at high confining stresses was explained by limited slip (most of the accommodation being done by bulk plastic deformation, a feature we did not introduce in our model, see above), and the lack of wear of smooth
surfaces is certainly related to an extremely small wear rate combined with a limited sliding distance. Hence, for what concerns gouge production, our model is much closer to the case of “rough” faults, despite it being nominally flat. This is explained by the presence of grain-scale irregularities which provide some interlocking in the initial state of the fault (Figure 6). A more realistic initial roughness would certainly improve the predictions of the model. To close this part of the discussion, it is useful to recall that wear rate estimates from field observations on real faults are generally much larger than those obtained in the lab: (Robertson, 1982) and (Scholz, 1987) reported wear rates (defined as the ratio between the thickness of the gouge layer and the slip distance) of 1,000 to 100,000 µm/m; (Katz et al., 2003) reported 1–20 mm of gouge thickness for sliding distances of 25–125 cm, leading to an average of ~10,000 µm/m, in satisfactory agreement with our study; and (Heesakkers et al., 2011) reported typical wear rates of 40,000 to 200,000 µm/m. These large values, which are closer to our numerical estimates, were possibly attributed to dynamic pulverization during one single M2.2 earthquake.

7.4. Spatial Variability on Faults

In the initial state of the system submitted to isotropic confining stress, prior to deviatoric loading, the fault can be considered as homogeneous. Some minor irregularities in the stress field are already present (Figure 11a), but the topography is nominally flat and the statistics of the discrete regions defining both surfaces (i.e., grain shapes, sizes and orientations, intact properties of the discretized rock, etc.) are the same at any location on the fault. We observe however that at the end of the first deviatoric elastic loading, a strong spatial variability of the stress field along the fault appears (Figure 11b). At that point, no major slip has occurred yet, but a large number of small events (which can be noticed as small stick-slip patterns at the end of the elastic part in Figure 11e and as sudden damage increases in the same time period in Figure 15) already took place at various locations of the fault. These events are local (i.e., they only rupture a small portion of the fault) and of small intensity (e.g., Figure 12c). This indicates that, from the very beginning of the simulation, some heterogeneities in the fault strength already exist: weak spots of several hundreds of µm yield and trigger local drops of the stress ratio σ<sub>nt</sub>/σ<sub>nn</sub> (as low as 0.25), thereby increasing the load on nearby strong spots of seemingly smaller spatial extensions (with endure stress ratios as high as 0.55). Variations in strength logically lead to variations in stress. At that scale of observation, the assumption of an Amontons-Coulomb friction law with a unique static friction coefficient on the whole fault proves wrong even before the first major slip event. The fact that S-faults have systematically larger first stress peaks than M-faults seems to indicate that such weak spots are more numerous on the latter, and that their occurrence is therefore more likely on larger faults. In contrast, some local regions on faults play the mechanical role of asperities (i.e., locus of shear stress concentrations) despite the absence of roughness. This certainly questions the usual concept of asperity, seen as a local geometrical feature belonging to the real-contact area, the latter being only a very small fraction of the whole interface. While this concept has been experimentally validated in the lab in the case of bare surfaces (Dieterich & Kilgore, 1994; Dillavou & Rubinstein, 2018), it is actually seldom used in daily industrial tribology, where third body accumulation quickly closes the gap between the surfaces and increases the proportion of real contact area in the interface to almost 1. We confirm this observation with the proposed model, but demonstrate that the concept of asperity defined as a localized stress concentration remains relevant in that case. The question of the lifespan of such an asperity across several seismic cycles remains open, though, and will be at the core of future works. It was indeed demonstrated based on thermal measurements in Aubry et al. (2018) that this lifespan is larger than the duration of a single large slip event, and this seems to be confirmed by a preliminary analysis of the present results (see e.g., the persistence of some stress patterns in Movie S2).

As the seismic cycles further accumulate, a large number of sliding events occur. Most of them are small or
intermediate (Figures 12b and 12c), and lead to local stress drops in the ruptured area and to stress increase at its
boundaries; And some of them are large (Figure 12a) and concern the whole fault, with sliding distances some-
times reaching a few tens of µm. Our results tend to indicate that the spatial distribution of the stress drop during
such events is also extremely variable (Figure 12a), and that they do not succeed in homogenizing the traction on
the fault. A quick comparison between Figures 11b and 11c indeed shows that the stress patterns before and after
the first major slip are similar in structure, at least to the first order. This observation is in good agreement with that
of Radiguet et al. (2013) who demonstrated numerically that (a) each rupture is affected by the state of stress of
the interface left by the previous ones, and that (b) the propagation of a single rupture front is not sufficient to
cancel an existing local stress concentration. The fault seems to rapidly reach a state of established disorder, with-
out any leveling of this disorder by further major events. Figure 11d however indicates that the spatial structure
of this disorder is prone to evolve in time, with an apparent increase of the correlation distance of the stress field as seismic cycles accumulate. Since the grain size is constant (no comminution is introduced in the code) we can speculate that the spatial scale driving this correlation distance may be the thickness of the gouge layer, which also increases in time. This hypothesis, however, would deserve to be validated by a more quantitative analysis in future works. A second scale of variation of the stress field is also visible on Figure 11d, corresponding to large wavelength variations of the shear stress associated with finite-strain effects at the boundaries of the system. Even in such an idealized system, we can therefore observe several nested scales of heterogeneities. In a recent paper (Mia et al., 2022) showed that off-fault plasticity may lead to partial ruptures as well as temporal clustering of seismic events in an initially homogeneous fault domain, through spatial redistribution of stress. This result is in line with the present findings. Additional complexity should also be expected from roughness, from heterogeneity in the rock properties or from interaction between nearby fault segments (Romanet et al., 2018). As proposed in McLaskey (2019), the common representation of a fault based on localized shear stress concentrations (asperities, strong fault patches, etc.) surrounded by a homogeneous ocean of lower stress is certainly oversimplified: in each weak or strong region, there might exist some weaker and stronger subregions, and so-on across the scales. Since rupture nucleation is often found to occur close to the boundaries between regions of different sliding behavior (e.g., at stress concentrations between locked and creeping regions (Chen & Lapusta, 2009)) we may wonder about the nature of the nucleation process when such boundaries are present literally everywhere on the fault. This may partly explain, for example, the wide range of nucleation lengths measured in McLaskey (2019) under nominally similar conditions but in the presence of an established gouge layer. As proposed in Schär et al. (2021), nucleation may occur by critical growth of a slip patch under sufficiently homogeneous conditions, but dynamic sliding in heterogeneous faults may also nucleate by progressive coalescence of countless yielding events at various scales, in a stepwise manner, following a “cascade up” logic (Noda et al., 2013).

## 8. Conclusion and Perspectives

The numerical results presented in this study demonstrate that a model combining a continuum-based and a discrete representation of rock can mimic the sliding behaviors observed in the laboratory earthquake experiments. These results are achieved based on a simple calibration of the bulk properties of the rock, without the need of ad-hoc weakening friction laws. Among the physical phenomena which are well-reproduced, one may cite the development of a large number of seismic cycles involving unstable slip events of various sizes and intensities, the progressive damaging of rock surfaces, and the production of a granular gouge layer in the interface. A remarkable complexity arises from this rather simple system, because of the limited randomness in the geometry of the discrete regions. The fields of normal stress, tangential stress, and effective friction in the vicinity of the sliding interface exhibit a considerable variability in space and time, with the spontaneous appearance of locked asperities based solely on granular phenomena, without the need of an associated geometrical feature. A quantitative analysis of the damage and wear of the interface reveal interesting trends: under the assumptions of the present model, we observe that the rate of damage is proportional to the accumulated strain in the loading system, that is, to the mechanical energy that is successively stored and released in the sample. Conversely, the thickness of the active gouge layer (participating to shear accommodation) increases at a different rate, and is proportional to the cumulated sliding distance. Both of these thicknesses increase with the confining stress. Quantitative estimates are broadly consistent with values inferred from field and lab studies.

The simulations depicted in this paper produced a large amount of data, of which only a small proportion was described above. These data include (a) the detailed motion and damage of each grain composing the gouge layer, (b) the detailed history of the kinetic and strain energies of the continuous regions of the sample, (c) the high-frequency (in space and time) record of the displacement and velocity fields of the surfaces bounding the discrete region of the sample, and (d) the same type of fields recordings for the full stress tensor. Future works will therefore be dedicated to exploit this data set in order to improve our understanding of the causes and consequences of the stress and strength heterogeneities we pointed out. For that purpose, we plan to explore the statistics of these heterogeneities: statistical moments and distributions, autocorrelation structures in space and time, etc. It may make it possible to evaluate what drives the characteristic scales of the strength and stress heterogeneities (Figure 11), and to determine to what extent these phenomena are modified by the progressive thickening of the gouge layer. The spatial distribution of damage in the fault zone can also be expressed in terms of roughness parameters (i.e., topography of the line separating intact and damaged rock), and this will be explored as well. Another angle of observation of these results is that of the sliding events. As we showed, the largest ones can be
detected from the stress drop they induce at the boundaries of the model (Figure 9), but most of them are way too small for such a detection (Figure 12c). We will therefore resort to the kinematic fields to identify all individual events of each simulation and characterize them quantitatively (hypocenter, spatial extension, duration, sliding distance, sliding velocity, complete energy budget, etc.). From this catalog we will explore the scaling relations between these events, and their possible clustering in foreshocks-mainshock-aftershocks causal families. We will also investigate their interactions with the heterogeneities in stress and strength, with particular focus on the way each event starts and stops. Finally, additional simulations may be considered in order to render the model closer to its experimental source of inspiration. This may include a more realistic roughness of the initial surfaces of the two half-samples, the introduction of heat production and diffusion (and its consequences on weakening and thermal cracking, including e.g., thermal pressurization), or the introduction of plasticity or poro-elasticity in the constitutive law of the bulk rock.

## Data Availability Statement

All simulations were performed with the open-source software MELODY version 3.94 (https://doi.org/10.5281/zenodo.4305614) developed by the first author and described in Mollon (2018). The code is available for downloading at https://github.com/GMollon/MELODY2D. The data generated is available at https://zenodo.org/record/7428795#.Y5dS44SZOF4.

## Acknowledgments

This work started during a sabbatical of the first author at the Laboratoire de Geologie de l'Ecole Normale Supérieure, in Paris, France, made possible by the CNRS. Both institutions are gratefully acknowledged.

## References

Aharonov, E., & Sparks, D. (2004). Stick-slip motion in simulated granular layers. *Journal of Geophysical Research*, 109(B9), B09306. https://doi.org/10.1029/2003jb002597

Archard, J. F. (1953). Contact and rubbing of flat surfaces. *Journal of Applied Physics*, 24(8), 981-988. https://doi.org/10.1063/1.1721448

Asadi, M. S., Rasouli, V., & Barla, G. (2012). A bonded particle model simulation of shear strength and asperity degradation for rough rock fractures. *Rock Mechanics and Rock Engineering*, 45, 649-675. https://doi.org/10.1007/s00603-012-0231-4

Aubry, J., Passelègue, F. X., Deldicque, D., Girault, F., Marty, S., Lahfid, A., et al. (2018). Frictional heating processes and energy budget during laboratory earthquakes. *Geophysical Research Letters*, 45(22), 12274-12282. https://doi.org/10.1029/2018gl079263

Aubry, J., Passelègue, F. X., Escartin, J., Gasc, J., Deldicque, D., & Schubnel, A. (2020). Fault stability across the seismogenic zone. *Journal of Geophysical Research: Solid Earth*, 125(8), e2020JB019670. https://doi.org/10.1029/2020jb019670

Bayart, E., Svetlizky, I., & Fineberg, J. (2018). Rupture dynamics of heterogeneous frictional interfaces. *Journal of Geophysical Research: Solid Earth*, 123(5), 3828-3848. https://doi.org/10.1002/2018jb015509

Belytschko, T., Lu, Y. Y., & Gu, L. (1994). Element-free Galerkin methods. *International Journal for Numerical Methods in Engineering*, 37(2), 229-256. https://doi.org/10.1002/nme.1620370205

Berthier, Y. (1990). Experimental evidence for friction and wear modelling. *Wear*, 139(1), 77-92. https://doi.org/10.1016/0043-1648(90)90210-2

Bolton, D. C., Shreedharan, S., Rivière, J., & Marone, C. (2020). Acoustic energy release during the laboratory seismic cycle: Insights on laboratory earthquake precursors and prediction. *Journal of Geophysical Research: Solid Earth*, 125(8), e2019JB018975. https://doi.org/10.1029/2019jb018975

Boneh, Y., Chang, J. C., Lockner, D. A., & Reches, Z. (2014). Evolution of wear and friction along experimental faults. *Pure and Applied Geophysics*, 171(11), 3125-3141. https://doi.org/10.1007/s00024-014-0801-3

Bowden, F. P., & Tabor, D. (1950). *The friction and lubrication of solids*. Oxford University Press.

Cao, H.-P., Renouf, M., Dubois, F., & Berthier, Y. (2011). Coupling continuous and discontinuous descriptions to model first body deformation in third body flows. *Journal of Tribology*, 133(4). https://doi.org/10.1115/1.4004881

Casas, N., Mollon, G., & Daouadji, A. (2022). DEM analyses of cemented granular fault gouges at the onset of seismic sliding: Peak strength, development of shear zones and kinematics. *Pure and Applied Geophysics*, 179(2), 679-707. https://doi.org/10.1007/s00024-021-02934-5

Chang, J. C., Lockner, D. A., & Reches, Z. (2012). Rapid acceleration leads to rapid weakening in earthquake-like laboratory experiments. *Science*, 338(6103), 101-105. https://doi.org/10.1126/science.1221195

Chen, T., & Lapusta, N. (2009). Scaling of small repeating earthquakes explained by interaction of seismic and aseismic slip in a rate and state fault model. *Journal of Geophysical Research*, 114(B1), B01311. https://doi.org/10.1029/2008jb005749

Cho, N., Martin, C. D., & Sego, D. C. (2008). Development of a shear zone in brittle rock subjected to direct shear. *International Journal of Rock Mechanics and Mining Sciences*, 45(8), 1335-1346. https://doi.org/10.1016/j.ijrmms.2008.01.019

Da Cruz, F., Emam, S., Prochnow, M., Roux, J. N., & Chevoir, F. (2005). Rheophysics of dense granular materials: Discrete simulation of plane shear flows. *Physical Review E*, 72(2), 1-17. https://doi.org/10.1103/physreve.72.021309

Dieterich, J. H., & Kilgore, B. D. (1994). Direct observation of frictional contacts: New insights for state-dependent properties. *Pure and Applied Geophysics*, 143(1-3), 283-302. https://doi.org/10.1007/bf00874332

Dillavou, S., & Rubinstein, S. M. (2018). Nonmonotonic aging and memory in a frictional interface. *Physical Review Letters*, 120(22), 224101. https://doi.org/10.1103/physrevlett.120.224101

Di Toro, G. D., Hirose, T., Nielsen, S., Pennacchioni, G., & Shimamoto, T. (2006). Natural and experimental evidence of melt lubrication of faults during earthquakes. *Science*, 311(5761), 647-649. https://doi.org/10.1126/science.1121012

Dorostkar, O., Guyer, R. A., Johnson, P. A., Marone, C., & Carmeliet, J. (2017). On the micromechanics of slip events in sheared, fluid-saturated fault gouge. *Geophysical Research Letters*, 44(12), 6101-6108. https://doi.org/10.1002/2017gl073768

Falk, M. L., Needleman, A., & Rice, J. R. (2001). A critical evaluation of cohesive zone models of dynamic fracture. In van der E. Giessen, S. Forest, & L. Kubin (Eds.), *Proceedings of the 5th European mechanics of materials conference on scale transitions from atomistics to continuum plasticity, March 5-8, 2001, Delft, The Netherlands, European mechanics of materials conference* (Vol. 11, pp. 43-50). Special Issue. Journal de Physique IV, Proceedings.
Feeny, B., Guran, A., Hinrichs, N., & Popp, K. (1998). A historical review on dry friction and stick-slip phenomena. *Applies Mechanics Reviews*, 51(5), 321–341. https://doi.org/10.1115/1.3099008

Ferrer, C., Salas, F., Pascual, M., & Orozco, J. (2010). Discrete acoustic emission waves during stick-slip friction between steel samples. *Tribology International*, 43(1–2), 1–6. https://doi.org/10.1016/j.triboint.2009.02.009

Fillot, N., Iordanoff, I., & Berthier, Y. (2007). Wear modeling and the third body concept. *Wear*, 262(7–8), 949–957. https://doi.org/10.1016/j.wear.2006.10.011

Fredrich, J. T., Evans, B., & Wong, T.-F. (1989). Micromechanics of the brittle to plastic transition in Carrara marble. *Journal of Geophysical Research*, 94(B4), 4129–4145. https://doi.org/10.1029/jb094ib04p04129

Gao, K., Euser, B. J., Rougier, E., Guyer, R. A., Lei, Z., Knight, E. E., et al. (2018). Modeling of stick-slip behavior in sheared granular fault gouge using the combined finite-discrete element method. *Journal of Geophysical Research: Solid Earth*, 123(7), 5774–5792. https://doi.org/10.1029/2018jb015668

Godet, M. (1984). The third-body approach: A mechanical view of wear. *Wear*, 100(1–3), 437–452. https://doi.org/10.1016/0043-1648(84)90025-5

Goebel, T. H. W., Schorlemmer, D., Becker, T. W., Dresen, G., & Sammis, C. G. (2013). Acoustic emissions document stress changes over many seismic cycles in stick-slip experiments. *Geophysical Research Letters*, 40(10), 2049–2054. https://doi.org/10.1002/grl.50507

Goldsby, D. L., & Tullis, T. E. (2011). Flash heating leads to low frictional strength of crustal rocks at earthquake slip rates. *Science*, 334(6053), 216–218. https://doi.org/10.1126/science.1207902

Guo, Y., & Morgan, J. K. (2004). Influence of normal stress and grain shape on granular friction: Results of discrete element simulations. *Journal of Geophysical Research*, 109(B12), B12305. https://doi.org/10.1029/2004jb003044

Guo, Y., & Morgan, J. K. (2007). Fault gouge evolution and its dependence on normal stress and rock strength—Results of discrete element simulations: Gouge zone properties. *Journal of Geophysical Research*, 112(B10), B10403. https://doi.org/10.1029/2006jb004524

Hazzard, J. F., & Young, R. P. (2000). Simulating acoustic emissions in bonded-particle models of rock. *International Journal of Rock Mechanics and Mining Sciences*, 37(5), 867–872. https://doi.org/10.1016/s1365-1609(00)00017-4

Heesakkers, V., Murphy, S., & Reches, Z. (2011). Earthquake rupture at focal depth, Part I: Structure and rupture of the Pretorius Fault, TauTona mine, South Africa. *Pure and Applied Geophysics*, 168(12), 2395–2425. https://doi.org/10.1007/s00024-011-0354-7

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*, 110(5), 1–18. https://doi.org/10.1029/2004jb003207

Houdoux, D., Amon, A., Marsan, D., Weiss, J., & Crassous, J. (2021). Micro-slips in an experimental granular shear band replicate the spatiotemporal characteristics of natural earthquakes. *Nature Communications Earth & Environment*, 2(1), 90. https://doi.org/10.1038/s43247-021-00147-1

Hsia, F.-C., Elam, F. M., Bonn, D., Weber, B., & Franklin, S. E. (2020). Wear particle dynamics drive the difference between repeated and non-repeated reciprocated sliding. *Tribology International*, 202, 105983. https://doi.org/10.1016/j.triboint.2019.105983

Iordanoff, I., Fillot, N., & Berthier, Y. (2005). Numerical study of a thin layer of cohesive particles under plane shearing. *Powder Technology*, 159(1), 46–54. https://doi.org/10.1016/j.powtec.2005.05.053

Kabir, M. A., Jasti, V. K., Higgs, C. F., & Lovell, M. R. (2008). An evaluation of the explicit finite-element method approach for modelling dense flows of discrete grains in a Couette shear cell. *Proceedings of the Institution of Mechanical Engineers—Part J: Journal of Engineering Tribology*, 222(6), 715–723. https://doi.org/10.1243/13506501jet388

Kandula, N., Cordonnier, B., Boller, E., Weiss, J., Dysthe, D. K., & Renard, F. (1999). Dynamics of microscale precursors during brittle compressive failure in Carrara marble. *Journal of Geophysical Research: Solid Earth*, 124(6), 6121–6139. https://doi.org/10.1029/2019jb017381

Katz, O., Reches, Z., & Baer, G. (2003). Faults and their associated host rock deformation: Part I. Structure of small faults in a quartz-syenite body, southern Israel. *Journal of Structural Geology*, 25(10), 1675–1689. https://doi.org/10.1016/s0191-8141(03)00011-7

Kim, H., & Buttlar, W. G. (2009). Discrete fracture modeling of asphalt concrete. *International Journal of Solids and Structures*, 46(13), 2593–2604. https://doi.org/10.1016/j.ijsolstr.2009.02.006

Latour, S., Schubnel, A., Nielsen, S., Madariaga, R., & Vinciguerra, S. (2013). Characterization of nucleation during laboratory earthquakes. *Geophysical Research Letters*, 40(19), 5064–5069. https://doi.org/10.1002/grl.50974

Lherminier, S., Planet, R., dit Vehel, V. L., Simon, G., Vanel, L., Máløy, K., & Ramos, O. (2019). Continuously sheared granular matter reproduces in detail seismicity laws. *Physical Review Letters*, 122(21), 218501. https://doi.org/10.1103/physrevlett.122.218501

Li, C.-Q., & Zhou, X.-P. (2021). Laboratory earthquake prediction of granite. *Tribology International*, 160, 107003. https://doi.org/10.1016/j.triboint.2021.107003

Mair, K., Frye, K. M., & Marone, C. (2002). Influence of grain characteristics on the friction of granular shear zones. *Journal of Geophysical Research*, 107(B10), ECV4-1-ECV4-9. https://doi.org/10.1029/2001jb000516

Mair, K., & Marone, C. (1999). Friction of simulated fault gouge for a wide range of velocities and normal stresses. *Journal of Geophysical Research*, 104(B12), 28899–28914. https://doi.org/10.1029/1999jb900279

Marty, S., Passelègue, F. X., Aubry, J., Bhat, H. S., Schubnel, A., & Madariaga, R. (2019). Origin of high-frequency radiation during laboratory earthquakes. *Geophysical Research Letters*, 46(7), 3755–3763. https://doi.org/10.1029/2018gl080519

McLaskey, G. C. (2019). Earthquake initiation from laboratory observations and implications for foreshocks. *Journal of Geophysical Research: Solid Earth*, 124(12), 12882–12904. https://doi.org/10.1029/2019jb018363

Mia, M. S., Abdelmeguid, M., & Elbanna, A. E. (2022). Spatio-temporal clustering of seismicity enabled by off-fault plasticity. *Geophysical Research Letters*, 49(8), e2021GL097601. https://doi.org/10.1029/2021gl097601

Mollon, G. (2015). A numerical framework for discrete modelling of friction and wear using Voronoi polyhedrons. *Tribology International*, 90, 343–355. https://doi.org/10.1016/j.triboint.2015.04.011

Mollon, G. (2016). A multibody meshfree strategy for the simulation of highly deformable granular materials. *International Journal for Numerical Methods in Engineering*, 108(12), 1477–1497. https://doi.org/10.1002/nme.5258

Mollon, G. (2018). A unified numerical framework for rigid and compliant granular materials. *Computational Particle Mechanics*, 5(4), 517–527. https://doi.org/10.1007/s40571-018-0187-6

Mollon, G. (2019). Solid flow regimes within dry sliding contacts. *Tribology Letters*, 67(4), 120. https://doi.org/10.1007/s11249-019-1233-0

Mollon, G., Aubry, J., & Schubnel, A. (2021). Simulating melting in seismic fault gouge. *Journal of Geophysical Research: Solid Earth*, 126(6), e2020JB021485. https://doi.org/10.1029/2020jb021485

Morgan, J. K. (1999). Numerical simulations of granular shear zones using the distinct element method: 2. Effects of particle size distribution and interparticle friction on mechanical behavior. *Journal of Geophysical Research*, 104(B2), 2721–2732. https://doi.org/10.1029/1998jb900055

Morgan, J. K., & Boettcher, M. S. (1999). Numerical simulations of granular shear zones using the distinct element method: 1. Shear zone kinematics and the micromechanics of localization. *Journal of Geophysical Research*, 104(B2), 2703–2719. https://doi.org/10.1029/1998jb900056
Nielsen, S., Spagnuolo, E., Smith, S. A. F., Violay, M., Toro, G. D., & Bistacchi, A. (2016). Scaling in natural and laboratory earthquakes. *Geophysical Research Letters*, 43(4), 1504–1510. https://doi.org/10.1002/2015gl067490

Noda, H., Nakatani, M., & Hori, T. (2013). Large nucleation before large earthquakes is sometimes skipped due to cascade-up-Implications from a rate and state simulation of faults with hierarchical asperities. *Journal of Geophysical Research: Solid Earth*, 118(6), 2924–2952. https://doi.org/10.1002/jgrb.50211

Okubo, K., Bhat, H. S., Rougier, E., Marty, S., Schubnel, A., Lei, Z., et al. (2019). Dynamics, radiation, and overall energy budget of earthquake rupture with Coseismic off-fault damage. *Journal of Geophysical Research: Solid Earth*, 124(11), 11771–11801. https://doi.org/10.1029/2019jb017304

Park, B., & Min, K.-B. (2015). Bonded-particle discrete element modeling of mechanical behavior of transversely isotropic rock. *International Journal of Rock Mechanics and Mining Sciences*, 76, 243–255. https://doi.org/10.1016/j.ijrmms.2015.03.014

Passelègue, F. X., Schubnel, A., Nielsen, S., Bhat, H. S., Deldicque, D., & Madariaga, R. (2016). Dynamic rupture processes inferred from laboratory microearthquakes. *Journal of Geophysical Research: Solid Earth*, 121(6), 4343–4365. https://doi.org/10.1002/2015jb012694

Passelègue, F. X., Spagnuolo, E., Violay, M., Nielsen, S., Di Toro, G., & Schubnel, A. (2016). Frictional evolution, acoustic emissions activity, and off-fault damage in simulated faults sheared at seismic slip rates. *Journal of geophysical Research: Solid Earth*, 121(10), 7490–7513. https://doi.org/10.1002/2016jb012988

Potyondy, D. O. (2007). Simulating stress corrosion with a bonded-particle model for rock. *International Journal of Rock Mechanics and Mining Sciences*, 44(5), 677–691. https://doi.org/10.1016/j.ijrmms.2006.10.002

Potyondy, D. O., & Cundall, P. A. (2004). A bonded-particle model for rock. *International Journal of Rock Mechanics and Mining Sciences*, 41(8), 1329–1364. https://doi.org/10.1016/j.ijrmms.2004.09.011

Quacquarelli, A. (2021). Modélisation multi-échelle de l'usure d'un outil diamanté pour l'optimisation microstructurale PhD Thesis. INSA Lyon.

Radiguet, M., Kammer, D. S., Gillet, P., & Molinari, J.-F. (2013). Survival of heterogeneous stress distributions created by precursory slip at frictional interfaces. *Physical Review Letters*, 111(16), 164302. https://doi.org/10.1103/physrevlett.111.164302

Rajagopal, P., Drozdz, M., Skelton, E. A., Lowe, M. J., & Craster, R. V. (2012). On the use of absorbing layers to simulate the propagation of elastic waves in unbounded isotropic media using commercially available Finite Element packages. *NDT & E International*, 51, 30–40. https://doi.org/10.1016/j.ndteint.2012.04.001

Reches, Z., & Lockner, D. A. (2010). Fault weakening and earthquake instability by powder lubrication. *Nature*, 467(7314), 452–455. https://doi.org/10.1038/nature09348

Renouf, M., Massi, F., Fillot, N., & Saulot, A. (2011). Numerical tribology of a dry contact. *Tribology International*, 44(7–8), 834–844. https://doi.org/10.1016/j.triboint.2011.02.008

Robertson, E. C. (1982). Continuous formation of gouge and breccia during fault displacement. In *23rd U.S symposium on rock mechanics (USRMS)*.

Romanet, P., Bhat, H. S., Jolivet, R., & Madariaga, R. (2018). Fast and slow slip events emerge. *Due to Fault Geometrical Complexity*, 45(10), 4809–4819. https://doi.org/10.1029/2018gl077579

Rouet-Leduc, B., Hulbert, C., Lubbers, N., Barros, K., Humphreys, C. J., & Johnson, P. A. (2017). Machine learning predicts laboratory earthquakes. *Geophysical Research Letters*, 44(18), 9276–9282. https://doi.org/10.1002/2017gl074677

Rubino, V., Lapusta, N., Rosakis, A. J., Leprince, S., & Avouac, J. P. (2014). Static laboratory earthquake measurements with the digital image correlation method. *Experimental Mechanics*, 55(1), 77–94. https://doi.org/10.1007/s11340-014-9893-z

Rubino, V., Rosakis, A. J., & Lapusta, N. (2017). Understanding dynamic friction through spontaneously evolving laboratory earthquakes. *Nature Communications*, 8(1), 15991. https://doi.org/10.1038/ncomms15991

Saksala, T., & Jabareen, M. (2019). Numerical modeling of rock failure under dynamic loading with polygonal elements. *International Journal for Numerical and Analytical Methods in Geomechanics*, 43(12), 2056–2074. https://doi.org/10.1002/nag.2947

Schär, S., Albertini, G., & Kammer, D. S. (2021). Nucleation of frictional sliding by coalescence of microslip. *International Journal of Solids and Structures*, 225, 111059. https://doi.org/10.1016/j.ijsolstr.2021.111059

Scholz, C. (1987). Wear and gouge formation in brittle faulting. *Geology*, 15(6), 493–495. https://doi.org/10.1130/0091-7613(1987)15<493:wagfib>2.0.co;2

Sobolev, G. A., Ponomarev, A. V., Koltsov, A. V., & Smirnov, V. B. (1996). Simulation of triggered earthquakes in the laboratory, induced seismic events (pp. 345–355).

Sone, H., & Shimamoto, T. (2009). Frictional resistance of fault during accelerating and decelerating earthquake slip. *Nature Geoscience*, 2, 750–708.

Vakis, A., Yastrebov, V., Scheibert, J., Nicola, L., Dini, D., Minfray, C., et al. (2018). Modeling and simulation in tribology across scales: An overview. *Tribology International*, 125, 169–199. https://doi.org/10.1016/j.triboint.2018.02.005

Wang, W., & Scholz, C. H. (1994). Wear processes during frictional sliding of rock: A theoretical and experimental study. *Journal of Geophysical Research*, 99(B4), 6789. https://doi.org/10.1029/93jb02875

Wang, Z., Wang, T., Wu, S., & Hao, Y. (2021). Investigation of microcracking behaviors in brittle rock using polygonal grain-based distinct method. *International Journal for Numerical and Analytical Methods in Geomechanics*, 45(13), 1871–1899. https://doi.org/10.1002/nag.3246

Wanne, T. S., & Young, R. P. (2008). Bonded-particle modeling of thermally fractured granite. *International Journal of Rock Mechanics and Mining Sciences*, 45(5), 789–799. https://doi.org/10.1016/j.ijrmms.2007.09.004

Xia, K., Rosakis, A. J., & Kanamori, H. (2004). Laboratory earthquakes: The sub-Rayleigh-to-supershear rupture transition. *Science*, 303(5665), 1859–1861. https://doi.org/10.1126/science.1094022

Xu, S., Fukuyama, E., & Yamashita, F. (2019). Robust estimation of rupture properties at propagating front of laboratory earthquakes. *Journal of Geophysical Research: Solid Earth*, 124(1), 766–787. https://doi.org/10.1029/2018jb016797

Yoshizawa, H., & Israelachvili, J. (1993). Fundamental mechanisms of interfacial friction. 2. Stick-slip friction of spherical and chain molecules. *Journal of Physical Chemistry*, 97(43), 11300–11313. https://doi.org/10.1021/j100145a031

Zhao, Z., Jing, L., & Neretnieks, I. (2012). Particle mechanics model for the effects of shear on solute retardation coefficient in rock fractures. *International Journal of Rock Mechanics and Mining Sciences*, 52, 92–102. https://doi.org/10.1016/j.ijrmms.2012.03.001