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

HAL
open science

# Laboratory Earthquakes Simulations-Emergence, Structure, and Evolution of Fault Heterogeneity

Guilhem Mollon, Jerome Aubry, Alexandre Schubnel

► To cite this version:

Guilhem Mollon, Jerome Aubry, Alexandre Schubnel. Laboratory Earthquakes Simulations-Emergence, Structure, and Evolution of Fault Heterogeneity. Journal of Geophysical Research: Solid Earth, 2024, 129 (6), 10.1029/2023JB028626. http://dx.doi.org/10.1029/2023JB0286273

HAL Id: hal-04662373

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

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_2.png" /></figure>

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

## RESEARCH ARTICLE

10.1029/2023JB028626

**Key Points:**

* A coupled discrete-continuum model of laboratory earthquakes allows the spontaneous development of fault complexity along seismic cycles

* Fault heterogeneity is persistent in space and time, with a strong positive correlation between normal stress, damage, and gouge thickness

* Stress fields exhibit a millimetric periodicity and a self-affinity at smaller scales, in quantitative agreement with previous estimates

**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. (2024). Laboratory earthquakes simulations—emergence, structure, and evolution of fault heterogeneity. *Journal of Geophysical Research: Solid Earth*, 129, e2023JB028626. https://doi.org/10.1029/2023JB028626

Received 22 DEC 2023
Accepted 9 JUN 2024

**Author Contributions:**

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

**Data curation:** Guilhem Mollon

**Investigation:** Guilhem Mollon

**Methodology:** Guilhem Mollon

**Software:** Guilhem Mollon

**Writing – original draft:**

Guilhem Mollon

**Writing – review & editing:**

Jérôme Aubry, Alexandre Schubnel

© 2024. The Author(s). 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_3.png" /></figure>

## Laboratory Earthquakes Simulations—Emergence, Structure, and Evolution of Fault Heterogeneity

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

¹Université de Lyon, INSA Lyon, CNRS, LaMCoS, Villeurbanne, France, ²Université Savoie Mont Blanc, Université Grenoble Alpes, CNRS, IRD, University Gustave Eiffel, Le Bourget du Lac, France, ³Laboratoire de Géologie, École Normale Supérieure/CNRS UMR 8538, PSL University, Paris, France

**Abstract** Seismic faults are known to exhibit a high level of spatial and temporal complexity, and the causes and consequences of this complexity have been the topic of numerous research works in the past decade. In this paper, we investigate the origins and the structure of this complexity by considering a numerical model of laboratory earthquake experiment, where we introduce a fault with homogeneous mechanical properties but allow it to evolve spontaneously to its natural level of complexity. This is achieved by coupling the elastic deformability of the off-fault medium (and therefore allowing for heterogeneous stress fields to develop) and the discrete degradation and gouge formation at the fault plane (and therefore allowing for structural heterogeneity to develop). Numerical results show the development of persistent stress, damage, and gouge thickness heterogeneities, with a much larger variability in space than in time. Strong positive correlations are found between these quantities, which suggest a positive feedback between local normal stress and damage rate, only mildly mitigated by the mobility of the granular gouge in the interface. For a wide range of confining stresses, after a sufficient number of seismic cycles, the fault reaches a state of established disorder with a constant roughness, a certain amount of periodicity at the millimetric scale, and a power law decay of the Power Spectral Density at smaller spatial scales. The typical height-to-wavelength ratio of geometrical asperities and the correlations between stress and damage profiles are in good agreement with previous field or lab estimates.

**Plain Language Summary** Earthquakes occur on faults in the Earth curst, and these faults are known to be highly complex, both from a geometrical and structural point of view. This complexity is now introduced artificially in many numerical and theoretical models of faults in order to derive meaningful knowledge related to earthquake mechanics. There is little knowledge, however, about the way this complexity develops. In this paper, we employ a numerical model to observe the spontaneous development of this complexity, in order to understand its origin, its structure, and its evolution.

### 1. Introduction

Seismic faults are broadly planar weakness zones in the lithosphere, comprising a relatively thin fault core surrounded by a damage zone, which is itself embedded in the intact country rock. The fault core is the region which accommodates most of the relative motion between the two sides of the fault plane, and is composed of fault rocks (gouge, cataclasite, breccia) that develop due to comminution or other processes which completely destroy the fabric of the initial rock (Tesei et al., 2014). The damage zone, in turn, is composed of country rock submitted to fracturing and secondary faulting, and keeps a memory of its initial fabric. Earthquakes generally occur along such faults which, being natural features associated with brittle weaknesses in the lithosphere, are subject to a large amount of complexity. One of these complexities is the spatial and temporal heterogeneity in fault properties that can be inferred from seismological measurements before, during, and after earthquakes (Scholz, 2019), as well as from field assessment (Chester et al., 1993; Mitchell & Faulkner, 2009; Shipton & Cowie, 2001) and geodetic studies (Perrin et al., 2016). A recent review of such heterogeneities is proposed in (Ben-Zion & Dresen, 2022; Tesei et al., 2014), including hierarchical geometrical heterogeneities and damage zones, that affect and complexify rupture dynamics (Beroza & Mikumo, 1996; Gounon et al., 2022; Latour et al., 2013). An important concept is the notion of asperity, seen as a local characteristic of the fault which is prone to concentrate phenomena such as rupture nucleation or initiation (Kemeny & Hagaman, 1992; Sammis et al., 1999; Selvadurai & Glaser, 2017; Zielke et al., 2017). In some contexts, asperities are seen as purely geometrical features 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 difficult to reconcile with field observations of natural faults, which are generally filled with a certain amount of gouge or cataclasite (Agosta & Aydin, 2006) and seem to have a real contact area close to the apparent one. In some other contexts, asperities are essentially regions of larger stress, either because of past sliding events or because of fault topography (Schmittbuhl et al., 2006). Finally, some studies represent fault asperities as purely rheological features, that is, as regions where the friction law is different (e.g., stronger, or with a different rate-dependency) from their surroundings (Chen & Lapusta, 2009; Dublanchet et al., 2013; Luo & Ampuero, 2018). Frictional heterogeneities have been reported to be related to a large number of local variations, such as fault fabric (Collettini et al., 2009), normal stress and fault roughness (Harbord et al., 2017), or rate of gouge production from initially bare surfaces (Noël et al., 2023). In fact, it is often stated that asperities are both geometrical and rheological features (Sagy & Brodsky, 2009), since geometry influences stresses (Kim et al., 2004), which in turn influence rheology. The precise causalities between these different fault properties (geometry, structure, rheology) remain however unclear, which raises the question of the emergence of fault heterogeneity. Another difficulty is related to scales since asperities are reported both at the regional scale (tens of km, Ghosh et al., 2012) and at the very local scale (tens of µm, Aubry et al., 2018).

Experiments on artificial heterogeneous faults (Bedford et al., 2022) showed that the average frictional response of such faults is not just an integrated average of the local friction properties. As a good example, a glimpse into the complexity of rupture propagation in a homogeneous fault with a localized patch of granular gouge may be found in Rubino et al. (2022). A recent trend in earthquake science is to investigate the consequences of fault complexity on phenomena relevant to earthquake mechanics, such as rupture nucleation, propagation, or arrest. Investigated heterogeneities cover properties alternation (Luo & Ampuero, 2018), roughness (Bruhat et al., 2020; Cattania & Segall, 2021; Tal et al., 2018), multiple segments (Romanet et al., 2018), heterogeneous stress level (Ripperger et al., 2007), weakening length (Albertini et al., 2021; Lebihain et al., 2021), off-fault inelasticity (Mia et al., 2022), variability in Rate-and-State Friction response (Chen & Lapusta, 2009; Dublanchet et al., 2013), etc. Many such studies rely on the artificial generation of random heterogeneous properties of the fault, to be introduced in numerical codes or in theoretical frameworks in order to derive meaningful conclusions. However, there is little knowledge about what should be the adequate properties of such random fields, for example, in terms of Probability Density Functions (PDF) or of autocorrelation functions of the spatially-varying quantities. A good way to observe the emergence of complexity is to perform laboratory earthquakes in controlled conditions, starting from a homogeneous initial state, and to monitor the spontaneous development of seismic cycles (Aubry et al., 2020; Bayart et al., 2018; Bolton et al., 2020; Dresen et al., 2020; Goebel et al., 2013; Guérin-Marthe et al., 2023; Kandula et al., 2019; Leeman et al., 2016; Li & Zhou, 2021; Marty et al., 2019; McLaskey, 2019; Passelègue et al., 2016; Scuderi et al., 2017; Sobolev et al., 1996; Xu et al., 2019). This can be nicely complemented by numerical simulations, provided that they are able to reproduce the main physics at stake, which include rock elastodynamics, damage mechanics, and granular physics (Casas et al., 2022, 2023; Dorostkar et al., 2017; Guo & Morgan, 2007; Mair et al., 2002; Mollon et al., 2021; Morgan & Boettcher, 1999; Papachristos et al., 2023; Taboada & Renouf, 2023; Wang et al., 2019). In the present paper, we revisit numerical data originating from a set of simulations described in details in Mollon et al. (2023). These simulations take inspiration from triaxial tests on sawcut marble samples reported in Aubry et al. (2018), albeit at a reduced scale and limited to 2D plane-strain kinematics (Figure 1). We employ this model to investigate and quantify in more depth this emerging fault complexity and provide hypothesis regarding its causes and consequences. Since they are inspired by a lab test, our simulations disregard some types of heterogeneities that are inherent to natural faults (e.g., lithological or large-scale geometrical ones), and focus on local quantities such as stress fields, gouge thickness, and damage intensity. Section 2 presents the main features of the model, and Section 3 describes the spatial variability of normal and shear stresses along the fault, and subsequently of the corresponding effective friction. It also focuses on the time variability of these quantities. In Section 4, the focus is put on the heterogeneity associated with the surface damage and the progressive build-up of the granular gouge layer, while Section 5 explores correlation structures in the data set. Section 6 proposes a discussion on the obtained results and links them to the literature on small-scale fault complexity, while Section 7 concludes this work.

## 2. Methods

The proposed model is based on a direct coupling between continuum and discrete representations of rock. The sample dimensions are 10 × 20 mm² (Medium-size, or “M-fault” in the nomenclature of Mollon et al. (2023), with a fault zone ~18 mm long), and its two halves are essentially represented as elastic solids, using a Lagrangian
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Coupled continuum-discrete numerical model of laboratory earthquakes. (a) Full view of the model and of its boundary conditions; (b) Zoom on the upper contact area between the sample and the loading system; (c) Zoom on a typical section of the fault; (d–f). Effective friction time series for confining stresses of 45, 90, and 180 MPa respectively. Insets are zooms on the early stages of the simulations, indicating the (approximate) time for yield $t_y$ and the first major stress drop $t_p$.</figcaption></figure>

meshfree continuum approach (Mollon, 2016, 2018). However, in the close vicinity of the fault plane (i.e., in the first 200 µm behind the surface on each side), a discrete representation is adopted where rock is represented as a collection of cohesively bonded rigid conformal polygons (Mollon, 2015; Saksala & Jabareen, 2019). This representation allows to simulate the progressive damage of the surface (breakage of the bonds between the grains) and the production of a gouge layer (kinematic separation of grains and transition toward a granular layer). In addition, an elastic bonding is used between the continuum regions and the discrete ones in order to ensure mechanical continuity. The elastic parameters of the continuum regions are chosen to simulate the Carrara marble used in the experiments (i.e., a Young's modulus of 70 GPa, a Poisson's ratio of 0.29, and a resulting shear modulus of 17.1 GPa). The contact parameters between the grains composing the discrete region follow a Cohesive Zone Model (CZM, Kim & Buttlar, 2009) which was calibrated by comparison with experimental results on Carrara marble taken from the literature (Fredrich et al., 1989), as described in details in Mollon et al. (2023). This calibration allows to represent in a fairly consistent way the elasticity of the intact rock (thanks to the elastic intact bonds between the grains) and its localized or distributed damage and breakage (thanks to the
elastic limit of these bonds and their post-breakage cohesion). The grain size is set to ~10 µm, which is much smaller than the typical size for calcite grains (~150 µm) in order to discretize as finely as possible these regions. This grain size corresponds to the typical size of large gouge particles identified in the lab (Passelègue et al., 2016) or in the field (Sagy & Brodsky, 2009). Overall, the model contains ~83,000 grains, as well as ~115,000 degrees of freedom in the elastic regions. An explicit integration scheme is used, with a time step of 0.2 ns. In the intact state of the discrete region, bonds are established between any pair of contacting grains. More specifically, a bond is created for any node of a grain A in contact with a segment of a grain B, and is attributed an initial damage equal to 0 and a length equal to the sum of the half-lengths of the two segments around the considered node of grain A. Along simulation, damage of each bond is only allowed to increase gradually from 0 to 1 following the CZM law described in Mollon et al. (2023). The damage attributed to each grain (and used in the metrics we employ in this work) is therefore a weighted sum of the damages of all its nodes. The model parameters are summarized in Table S1 in Supporting Information S1.

The two half-samples are initially separated, and positioned between two deformable elastic blocks (Figure 1), which simulate the compliance of the loading system and allow a better dissipation of the elastic waves traveling in the sample. A confining stress $\sigma_3$ is first applied on the sides of the sample and on the loading system, which activates the contacts between the half-samples along the fault plane. This stage is pursued until mechanical equilibrium is reached. The boundary conditions of the loading system are then switched to displacement-driven: the lower boundary of the system is maintained fixed while the upper boundary is assigned a downwards velocity which leads to an average vertical strain rate of $2.5 \text{ s}^{-1}$ in the whole system. This loading rate is orders of magnitude faster than the experimental one, and is chosen as a trade-off between computational costs and system stick-slip dynamics (Mollon et al., 2023; Scott, 1996): it is chosen as high as possible, provided that it allows to clearly identify dynamic sliding events and locked periods. Three confining stresses of 45, 90, and 180 MPa are tested (labeled M45, M90, and M180 respectively), as in the target experiments. As shown in the curves of Figures 1d–1f, all simulations initiate with an elastic loading period, where shear stress on the fault plane increases steadily until a yield point (time $t_y$). The curve slope then starts to slowly decrease, and fluctuations in the shear stress appear (See insets in Figure 1.). At the end of this yield period, at time $t_p$, a first major stress drop occurs, with a peak effective friction close to 0.3 for all three confining stresses. From thereon, all three faults enter into a succession of a large number of seismic cycles, with successive episodes of fault locking and sliding. The yield point can only be roughly estimated because yield occurs gradually, but we extract $t_y \approx 0.13, 0.31$, and 0.39 ms for confining stresses of 45, 90, and 180 MPa respectively. In contrast, the times for the first major stress drops are well-defined as $t_p = 0.32, 0.58$, and 0.96 ms respectively. As illustrated in Mollon et al. (2023), the inter-event times, durations, average velocities, sliding distances, and stress drops of these events are extremely variable: the largest events rupture the whole 18 mm of the fault and reach both extremities, leading to stress drops of several tens of MPa and sliding distances of several tens of µm; while the smallest events are entirely localized in small patches inside the fault without reaching its extremities, and are not even detectable as stress drops in the loading system. This variability is made possible by the elastic compliance of the rock surrounding the fault, in contrast with more common Discrete Element Modeling (DEM) where rigid walls impose a perfectly synchronized sliding along the whole simulated fault zones (see a discussion in Papachristos et al., 2023).

No precise statistics were extracted yet from this catalog of events, and this work is kept for future studies since it is not at the core of the present paper. This stick-slip accommodation regime is pursued until a total sliding distance of ~500 µm is reached for each of the three simulations, which corresponds to a simulated time of 5 ms. During all this duration, and as hundreds of seismic cycles of various sizes occur, the contacting surfaces get more and more damaged, grains detach from the damaged rock, and a gouge layer is progressively formed. This may modify the frictional response of the fault by delocalizing slip accommodation, and indeed the first few major stress drops are more intense than the following ones because they occur at larger peak shear stresses, before a fully granular gouge layer is established. The initial thickness of the discrete region (200 µm on each surface) is sufficient to avoid wall-effects, which means that the “damage front” which progresses deeper and deeper along time does not reach the interfacial zone between the discrete and continuum regions. As shown and quantified in Mollon et al. (2023), the rate of thickening of this damage zone is proportional to the confining stress and its thickness eventually increases linearly with the cumulated strain in the sample. In contrast, the rate of thickening of the gouge layer participating to accommodation remains largely independent from the confining stress, and the gouge thickness eventually increases linearly with the cumulated sliding distance on the fault plane. These results are only in partial agreement with the Archard law (Archard, 1953):
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Normal $\sigma_{nn}$ and tangential $\sigma_{nt}$ stress heterogeneities on the upper half of the numerical fault (hanging wall) under a confining stress $\sigma_3 = 90$ MPa. Black dashed line indicates the first major sliding event, which can also be noticed in the friction history on the left.</figcaption></figure>

the proportionality of wear volume to both sliding distance and normal stress is not fully ensured, and the very notion of wear rate is rather ill-defined (Mollon et al., 2023). Interestingly, the Archard law is now often described as inaccurate and oversimplified in the tribological literature (Hu et al., 2022; Reichelt & Cappella, 2021). It is worth mentioning at that point that the simulated fault is mechanically homogeneous in its initial state. It is nominally flat, apart from a small-scale roughness related to the shapes of the surface grains, and the initially bonded granular packing has homogeneous statistics in terms of sizes, shapes, and orientations of the grains (Figure 1c). Yet, numerical results provided in Mollon et al. (2023) indicate that fault heterogeneity naturally emerges from this simulated system, under the form of localized slip events, patches of stress concentration, or localized damaged areas.

### 3. Variability in Stress and Friction

The continuum elastic character of the major part of the sample makes it possible to monitor the full stress tensor field during the entire duration of a simulation (Mollon et al., 2023). It is useful to project this stress tensor in the local frame of the fault plane $\{\vec{n}, \vec{t}\}$ (Figure 1a) and to consider the normal stress $\sigma_{nn}$ and the shear stress $\sigma_{nt}$ expressed in that frame. Such fields are plotted in the upper part of Figure 2, in the upper half-sample at the end of the simulation M90. Two levels of heterogeneity are observed. At the scale of the whole sample, both normal and shear stress are subject to large variations which are attributed to finite-strain effects (Mollon et al., 2023). More specifically, at that stage of the simulation the fault has slid by ~500 µm, while the horizontal motion of the upper and lower half samples were restricted by their frictional contact with the loading system (in accordance with the target experimental conditions), leading to this loss of macroscopic stress homogeneity.

At a smaller scale and in the vicinity of the fault plane, a different kind of heterogeneity is observed, which is directly related to phenomena occurring in the fault interface during the seismic cycles. The lower part of Figure 2
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Evolution of the shear stress field and of the damage distribution in a selected area of the numerical fault (marked as a white rectangle in the upper-left caption) under a confining stress $\sigma_3$ = 90 MPa. Large-scale stress heterogeneities in the sample (upper-left caption) are due to lateral kinematic restrictions at the boundaries of the sample, in accordance with the experiments where these boundaries were not lubricated, as commented in (Mollon et al. (2023)).</figcaption></figure>

(and Figure S1 in Supporting Information S1) allows to follow these heterogeneities in time, alongside the resulting average friction measured on the fault. Both normal and shear stress profiles were first filtered (gaussian filter with a sigma-value of 80 µm) in order to remove small scale noise related to the local contact conditions at the boundary between continuous and discrete regions. We observe that both normal and shear stresses develop large heterogeneities along the fault plane. These heterogeneities appear to grow progressively during the initial elastic loading, that is, before the first major stress drop, and to keep growing during the first few seismic cycles, up to $t = 1$ ms (especially $\sigma_{nn}$). When slip accumulates, stress patterns appear to develop and to remain somewhat persistent in time. The average normal stress on the fault is close to 115 MPa, but local values can range between 20 and 200 MPa. Likewise, the average shear stress is close to 28 MPa (with an average friction on the fault close to 0.24), but typically ranges between 5 and 60 MPa. Such large variations of the interfacial stresses are rather unexpected on such a small system with initially homogeneous mechanical properties.

This is illustrated in more details in Figure 3 (completed by Supplementary Movie S1), which focuses on a selected patch on the fault during the whole M90 simulation and provides both the field of damage in the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Illustration of the spatial shift between the normal and tangential stresses measured at the interface between the continuous and discrete regions ($\sigma_3$ = 90 MPa, $t$ = 5 ms, end of the simulation): (a) Normal stress concentration; (b) Shear stress concentration at the same location; (c) Orientations of the normal and shear stress patterns; (d) Cross-correlation function between normal and shear stress on the hanging wall (computed on the whole fault and during the entire duration of the simulation), with a clear maximum for a spatial lag of 280 µm.</figcaption></figure>

degradable region and the field of shear stress around it. Several interesting observations can be made. At $t$ = 0.2 ms, the fault is still in a state of elastic loading, and no evidence of damage is visible. During the yield period, at $t$ = 0.4 ms, first evidences of damage start to appear way before the first major stress drop (which occurs around $t$ = 0.58 ms in that case), with the appearance of a thin and discontinuous layer of damaged grains in the interface and some hints of localized crack initiations in the degradable region of the sample. Moderate shear stress concentrations can be observed on each side of the discrete region. At $t$ = 0.6 ms, a first major sliding event just took place, and the damaged layer is now continuous along the whole fault. We also notice that this layer is not perfectly straight (i.e., we can observe a surface at the gouge-host rock interface rougher than the initial bare surface), with some more intense cracking patterns initiating at several locations in the interface and progressing toward the bulk intact rock. These cracks seem to correlate well with the locations of the shear stress concentrations, which increase in amplitude (see for example the lower part of the observation window). A direct comparison of these damaging patterns with experimental ones is rendered difficult by the relatively large size of the discrete grains, which does not allow a sufficient resolution to derive conclusions on such phenomena as crack preferred orientation (Renard et al., 2018). This resolution is however sufficient to describe the main phenomena at stake, that is, progressive damaging and wear. As sliding progresses and as seismic cycles accumulate, the central layer of fully damaged grains gets thicker, which corresponds to the progressive build-up of a gouge layer. Meanwhile, a damage front progresses toward intact rock, but not in a homogeneous manner. At $t$ = 2.0 ms, for example, we can clearly discern two areas of larger damage (associated with an intense crack network) separated and surrounded by three patches of more intact rock. The shear stress patterns also gain in intensity, and sometimes fluctuate in time, but remain essentially stationary on the fault despite the large number of sliding events that take place between each caption. At $t$ = 5.0 ms, we observe that the damage front is getting close to the interface between the discrete and continuum regions, which indicates that simulations should be stopped in order to avoid wall-effects. Generally speaking, it seems that shear stress concentrations tend to coalesce and to be less numerous and more distant as the fault matures (see also Figure 2). This can be clearly seen by comparing the stress patterns at $t$ = 1.5 ms, where two areas of stress concentration are visible on each side of the fault, and $t$ = 5.0 ms, where only one large stress concentration is visible on each side of the fault in this area. It seems reasonable to call this area an asperity, since it appears to be permanently established (Figure 2) and to bear a much larger load than its surroundings.

It is interesting to observe that, because of the finite thickness of the degradable region (~400 µm, including the granular gouge layer), the shear stresses measured on the footwall and on the hanging wall (at the boundary between the continuous and discrete regions) are not equal (Figure 3, top-left). This is not related to a lack of mechanical equilibrium, but to the fact that they are measured at a certain distance from the central “mathematical” plane of the fault (on which continuity of stress fields is guaranteed). Hence, if the shear stress history on the foot wall was plotted in the same manner as in Figure 2, it would not coincide with that plotted for the hanging wall (although the integrals of these stresses on the entire fault are exactly balanced when the fault is locked). As a matter of fact, under such spatial fluctuations of the shear stress and in the presence of a relatively thick gouge layer, it is difficult to define a univocal value of the shear stress at a given location of the fault plane, and this concept is only valid in idealized models based on an infinitely thin fault zone (Wibberley et al., 2008). Figure 4 provides a way to solve this
issue. It shows, at the same location as in Figure 3, the fields of normal stress (Figure 4a) and of shear stress (Figure 4b), and seems to indicate that the normal stress is symmetrically distributed while the shear stress is shifted between the footwall and the hanging wall. This is confirmed quantitatively in Figure 4d, which shows the cross-correlation function between the normal and shear stress on the hanging wall. A clear peak is observed for a spatial shift of 280 µm, which indicates that $\sigma_{nn}$ and $\sigma_{nt}$ are strongly correlated in space (0.64) if this fault-parallel shift of the shear stress is accounted for. Histograms of Figure S2 in Supporting Information S1 provide additional evidence for this spatial correlation. Hence, in order to evaluate the local value of the fault effective friction, a relevant technique is to consider the sketch of Figure 4c. The normal stress at point O can be approximated by $\sigma_{nn,O} = (\sigma_{nn,M} + \sigma_{nn,N})/2$, and the shear stress at point O can be approximated by $\sigma_{nt,O} = (\sigma_{nt,P} + \sigma_{nt,Q})/2$. The effective friction $\mu_O$ at point O is therefore computed the following way:

$$\mu_O = (\sigma_{nt,P} + \sigma_{nt,Q})/(\sigma_{nn,M} + \sigma_{nn,N}) \quad (1)$$

where $\sigma_{nt,P}$ and $\sigma_{nt,Q}$ are the shear stresses measured at points P and Q (shifted by ±280 µm with respect to O, a value which is consistent for all our simulations), and $\sigma_{nn,M}$ and $\sigma_{nn,N}$ are the normal stresses measured at points M and N (aligned with O).

In order to visualize the variety of the stress states along the fault, Figure 5 shows bivariate histograms of the distribution of $\sigma_{nn}$ and $\sigma_{nt}$ as expressed in the previous paragraph. In Figures 5a–5c, dashed lines correspond to constant ratios of $\sigma_{nt}$ and $\sigma_{nn}$, that is, of constant effective friction, ranging between 0.1 and 0.5. These distributions confirm that normal and shear stresses are strongly correlated, since the bivariate histograms form very elongated shapes in all cases. These shapes do not align with the origin, which indicates that they do not correspond to a given friction coefficient. They do not align either with the theoretical loading path (white thick lines) that would be expected on the fault plane in a perfectly homogeneous case:

$$\sigma_{nn} = \frac{\sigma_1 + \sigma_3}{2} + \frac{\sigma_1 - \sigma_3}{2} \cos(2\theta) \quad (2)$$

$$\sigma_{nt} = \frac{\sigma_1 - \sigma_3}{2} \sin(2\theta) \quad (3)$$

where $\theta = 60^\circ$ is the angle between the fault plane and the minimum principal stress $\sigma_3$ (which is equal to 45, 90, or 180 MPa), and $\sigma_1$ is allowed to vary. An analysis of the motion of each local stress state in the $\sigma_{nn}-\sigma_{nt}$ frame over the seismic cycles requires the precise definition of a catalog of events and the clear identification of the stress state at failure, at each location and for each event. This work is kept for a future study, but for now, it is interesting to observe (Figure 5d) that the local effective friction is generally comprised between 0.1 and 0.3, which seems to be a point of fault “saturation” at all confining stresses. It is also interesting to point out that the relative variability in stress states seems to be larger at low confining stresses (M45, Figure 5a) but more limited under larger confinement (M180, Figure 5c).

A map of the effective friction on the fault along time is provided in Figure 6a for the case $\sigma_3 = 90$ MPa. It shows very clear patterns in space, with strong oscillations along the fault that remain generally persistent in time. When averaged over the post-peak period of the simulation (Figure 6b), these oscillations are found to remain between 0.15 and 0.3, and a periodicity clearly appears. This is confirmed by the average autocorrelation function of the effective friction (Figure 6c), which exhibits a very clear secondary peak for a distance of 0.94 mm. Similar results are obtained for the cases M45 and M180 (Figure S3 in Supporting Information S1), with respective spatial periods of 1.06 and 0.96 mm. It therefore appears that this periodicity is independent from the confining stress and that the effective friction period remains close to 1 mm in all simulations. An interesting point is that these friction patterns seem to “nucleate” during the pre-peak period in a much clearer manner than the stress patterns observed in Figure 2. These initial heterogeneities in the elastic regime could arise (a) from the small irregularities in the intact surfaces, related to the geometry of the grains, or (b) from the local patterns in the intact microstructure of the discrete region. Unfortunately, it is not possible with the data at hand to discriminate between these two potential sources of heterogeneity. We can be sure, however, that neither the continuous region not its connection to the discrete region is the source of these patterns, as its homogeneity is nearly perfect.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Stress state distribution on the faults; (a–c). 2D histograms showing the joint distributions of normal and tangential stresses on the fault at each location and time of the simulations (ignoring the elastic loading period before the first major stress drop). Dashed white lines indicate effective friction coefficient and thick solid lines indicate the theoretical loading path on the fault plane; (d) The three histograms superimposed, with additional black lines highlighting the major data structures.</figcaption></figure>

In order to get a better picture of the structure of the stress heterogeneities, it is useful to distinguish the stress
variabilities in time and in space. This is the purpose of Figure 7, which focuses on the field of normal stress on the
hanging wall of the fault (identical results were obtained for the foot wall). In Figure 7a, we focus on the vari-
ability in time: at each time increment, the mean value and the standard deviation of the normal stress σₙₙ are
computed, and are plotted for the three cases M45, M90, and M180. The three simulations follow similar trends,
despite some minor differences: during the elastic loading, we observe a steady increase of the space-averaged
normal stress, with a constant standard deviation. At the first major stress drop, both the average and the stan-
dard deviation decrease, which means that this event slightly homogenizes the stress state along the fault.
However, as the seismic cycles accumulate, we generally observe an increase in the standard deviation of σₙₙ. In
the cases M45 and M180 this spatial variability tends to stabilize after ~3 ms of simulation (i.e., ~300 µm of
sliding), while it keeps increasing more steadily until the end of the simulation for the case M90. We observe that,
despite the rather high level of spatial variability, the spatial average of σₙₙ is not submitted to large fluctuations in
time, and generally increases very slowly during all simulations.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Effective friction statistics for a confining stress of 90 MPa; (a) Space-time map of the effective friction along the fault. The black dashed line indicates the first major stress drop; (b) Time-averaged friction curve; (c) Autocorrelation function of the effective friction, (b) and (c) ignore the elastic loading period preceding the first major stress drops.</figcaption></figure>

Figure 7b (completed by Supplementary Movie S2) adopts the opposite perspective: for each location on the fault, a mean value of $\sigma_{nn}$ is computed for the time period running from the first major stress drop to the end of the simulation (marked in light color in Figure 7a). Likewise, a variability in time is computed through a standard deviation of $\sigma_{nn}$ on the same period. The obtained curves reveal that the time-averaged normal stress strongly

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Normal stress heterogeneities. (a) Evolution over time of the spatial average and standard deviation of $\sigma_{nn}$; (b) Variation along the fault of the time-average and standard deviation of $\sigma_{nn}$. Both (a) and (b) are compared to the general average and standard deviation computed on the whole simulation (excluding initial elastic loading part).</figcaption></figure>
fluctuates in space, and that the amplitude of these fluctuations is much larger than the variation of $\sigma_{nn}$ in time. In the case M90, for example, the peak-to-trough amplitude in space is close to ~60 MPa, while the typical standard deviation in time at a given location is of the order of ~15 MPa. Apparently, both space and time variabilities increase with the confining stress. Generally speaking, this figure draws the picture of a fault subjected to persistent patterns of normal stress concentration with limited fluctuation in time, despite the large number of sliding events of all sizes taking place in the interface. This confirms the visual impression derived from Figure 2.

4. Variability in Damage and Gouge Thickness

To investigate the origins of the stress variability described in the previous sections, it is instructive to focus on the phenomena that occur locally in the interface. As shown in Mollon et al. (2023) and illustrated in Figure 3, the fault plane is subjected to a progressive damage of the initially intact rock, and to the production of a granular gouge which progressively accumulates in a central layer. Despite the initial homogeneity and planarity of the rock surfaces, Figure 3 shows that we should expect a spatial variability regarding these phenomena. In order to quantify this, we adopt the procedure described in Figure 8. Figure 8a shows a selected area of the fault plane, where each grain (may it belong to intact rock, damaged rock, or gouge) is colored with respect to its level of relative damage. A relative damage equal to 0 means that all the bonds this grain had with initial neighbors are still intact, while a relative damage of 1 means that all these bonds were broken. We then interpolate damage profiles in the direction orthogonal to the fault plane every 20 µm (Figure 8b), and extract the transversal positions corresponding to relative damages of 0.99 and 0.4, which gives four profiles along the fault (two profiles on each side of the fault plane). These profiles are gaussian-filtered (sigma-value of 80 µm), and are then used to define two quantities at each location of the fault: (a) the gouge thickness $T_G$ (region comprised between the two 0.99 profiles), and (b) the damaged thickness $T_D$ (sum of the regions comprised between the 0.4 and 0.99 profiles on each side, Figure 8c).

These profiles evolve in time as seismic cycles accumulate, and the resulting space-time maps are provided in Figure 9 (Completed by Supplementary Movie S3) for a confining stress $\sigma_3 = 90$ MPa (Figure S4 in Supporting Information S1 provides the same data for the other cases). Since the damaging process is irreversible under the assumptions of the present model (no healing is implemented), $T_G$ and $T_D$ can only increase in time, and indeed this is what is observed in Figure 9. More specifically, the damaged thickness starts increasing early in the simulation, at a time which we can broadly associate with the yield point $t_y$, and well-before the first major sliding event. This is in good agreement with observations provided in Figure 3. Conversely, the gouge thickness starts increasing later, when approaching the first major event. When this first large stress drop occurs at $t_p$, important spatial heterogeneities start to develop, both for $T_G$ and $T_D$. As was observed for stresses (Figure 2) and friction (Figure 6), these heterogeneities are persistent in time despite the continuous shearing of the interface. A clear correlation appears between the spatial distribution of $T_G$ and $T_D$, although the regions of higher damage seem “sharper” (i.e., less spatially extended) than those of larger gouge thickness. This is certainly to be attributed to the larger mobility of the grains within the gouge layer, which allows them to spread along the fault and to smoothen the gouge thickness profile. This smoothing effect is not sufficient, however, to produce a perfectly homogeneous gouge layer. Similar observations can be derived from Figure S4 in Supporting Information S1, for the cases M45 and M180.

A more quantitative picture is drawn in Figure 10, which provides the evolution in time of key quantities. Figures 10a and 10b show the evolution in time of the average gouge and damage thicknesses for the three simulations (M45, M90, and M180), along with their standard deviations in space. As mentioned in Mollon et al. (2023), it shows that the damaging process follows different stages: after a short period of fully elastic loading without damage, this process strongly accelerates because of the local breakage of small asperities and interlocking grains, until a large jump is observed during the first major sliding event. The rate of damage and of gouge production then progressively decreases as a thicker gouge layer protects the surfaces, and seems to stabilize to a constant, low damage rate toward the end of the simulations (Boneh et al., 2014; Queener et al., 1965). Figure 10a indicates that the gouge layer is slightly larger when the confining stress is increased, but that its spatial variability tends to reach a constant level. The spatial standard deviation of gouge thickness (as materialized by the width of the colored patches around the mean curves) stabilizes close to 13 µm for all three confining pressures. Faults therefore seem to converge toward a constant state of disorder in terms of gouge thickness, despite a steady increase of the mean value of this thickness. Conversely, Figure 10b shows that the spatial variability in damage thickness increases both with the sliding distance and with the confining stress: from the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Determination of the distribution of the gouge and damage thicknesses $T_G$ and $T_D$. (a) Close-in view of a selected region of the fault with each grain colored according to its relative damage state; (b) Interpolation of transversal damage profiles at each location along the fault; (c) Thresholding of the damage distribution and definition of the gouge and damaged regions (unfiltered on this figure).</figcaption></figure>

first major stress drop to the end of the simulation, the standard deviation of $T_D$ increases from 9 to 15 µm for $\sigma_3$ = 45 MPa, and from 15 to 40 µm for $\sigma_3$ = 180 MPa. It therefore appears that the fault complexity in terms of damage thickness increases with slip, and that this increase is enhanced at larger confining stresses. Combined with the persistence in space of damage patterns observed in Figure 9, this result tends to point toward a positive feedback of local damage, that is, to the fact that areas with a larger damage at a given time are likely to exhibit a larger rate of damaging in later stages of the simulations.

Figure 10c shows the evolution of the fault quadratic roughness $R_q$ (root-mean-square of the heights), measured at the boundaries between the gouge layer and the damaged rock (i.e., at the profiles corresponding to a relative damage equal to 0.99, Figure 8). For all the confining stresses, it shows a strong and steady increase of the roughness during the yield period (i.e., between $t_y$ and $t_p$), in concomitance with the initiation of damage.
<figure><img src="image_12.png" /><figcaption>Figure 9. Space-time maps of the evolution of the gouge thickness $T_G$ and the damage thickness $T_D$ under a confining stress $\sigma_3 = 90$ MPa.</figcaption></figure>

However, this roughness peaks when the first major sliding event occurs, and then does not vary very much during the following numerous seismic cycles: it slowly decreases toward a constant value for M45 and M180, but steadily increases for M90. The numerical value of quadratic roughness seems to be hardly affected by the confining stress: it peaks between 8 and 9 µm for all cases, and stabilizes between 7 and 9 µm at the end of the simulations. Hence, in a similar way to the gouge thickness, faults appear to tend toward an established roughness, roughly independent from the sliding distance and the confining stress, with a quadratic roughness slightly lower than one grain diameter (equal to 10 µm under the assumptions of the present model). Fault profiles at early stages (i.e., just before the first major event) and in a mature state (at the end of the simulations) are provided in Figures 10d and 10e, and confirm visually the previous observations: The progression of the damage front toward intact rock is faster and more heterogeneous for large confining stresses, but this trend is not reproduced by the gouge layer profile because the larger grain mobility favors reorganization and therefore attenuates these heterogeneities.

## 5. Correlation Structures

Figures 2, 6, and 9 clearly demonstrated the existence of patterns in stress, friction, and damage along the fault. The relationship between these patterns, however, deserves a deeper analysis. In this section, we focus on three quantities: $\sigma_{nn}$, $T_G$, and $T_D$. As illustrated in Figures 7 and 10, these quantities all vary in time following some increasing trend. In order to decouple the spatial patterns from the temporal evolution, we define the quantities $\Delta\sigma_{nn}$, $\Delta T_G$, and $\Delta T_D$ in the following way:

$$ \Delta\sigma_{nn}(x,t) = \sigma_{nn}(x,t) - \overline{\sigma_{nn}(t)} \quad (4) $$

where $\overline{\sigma_{nn}(t)}$ is the spatial average of $\sigma_{nn}(x, t)$ at a given time $t$, and likewise for $\Delta T_G$ and $\Delta T_D$. The PDF of these three quantities are provided in Figure 11 for the simulations M45, M90, and M180. The deviation in normal stress $\Delta\sigma_{nn}$ is normalized by the applied confining stress $\sigma_3$. These distributions remain essentially symmetric, but are bounded and have a somewhat triangular shape, in particular for the distribution of $\Delta\sigma_{nn}$. It appears that a large confining stress tends to reduce the relative variability of the normal stress on the fault plane, but to increase the variability in damage thickness. The variability of gouge thickness seems relatively unaffected.

We can then analyze the possible correlations between these three quantities. The results are provided in Figure 12 for the case M90, and in Figure S5 in Supporting Information S1 for M45 and M180. Figure 12 shows the
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Evolution of the fault complexity. (a) Evolution of the average gouge thickness (spatial average ± one standard deviation); (b) Evolution of the average damage thickness (spatial average ± one standard deviation); (c) Evolution of the fault quadratic roughness $R_q$; (d) Fault complexity at an early stage (end of the elastic loading, prior to stress peak, at $t = t_p - 0.01$ ms); (e) Fault complexity at a late stage (end of the simulation). For (d) and (e), the vertical scale is dilated by an arbitrary factor 10 to ease visualization. The gouge layers are represented in light yellow, and the damaged layers in orange.</figcaption></figure>

statistical distribution of every fault state (i.e., each location at each time of the simulation, ignoring the elastic loading period before the first major stress drop) in the 3D space {$\Delta\sigma_{nn}$, $\Delta T_G$, $\Delta T_D$}. This distribution is represented by three 2D histograms in the three planes defined by each pair of quantities. They are all centered on 0 by construction. A zero-mean trivariate normal distribution is fitted on this empirical distribution, and

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Probability Density Functions (PDF, arbitrary units) of the spatial deviations in normal stress (normalized by the confining stress), in gouge thickness, and in damage thickness, for the three confining stresses. Elastic loading before the first major stress drop is ignored in these statistics.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Correlation structure of the spatial deviations $\Delta\sigma_{nn}$, $\Delta T_G$ and $\Delta T_D$ at each location and time of the numerical fault, illustrated by 2D histograms and isodensity ellipsoids of the best-fitting trivariate normal distribution projected in the plane of each pair of quantities, for a confining stress of 90 MPa. Each figure also indicates (i) the slope of the largest principal axis of this ellipsoid in each 2D plane and (ii) the coefficient of correlation between the quantities plotted in each plane. Elastic loading before the first major stress drop is ignored in these statistics.</figcaption></figure>

represented in the same planes. It is found to fit relatively well the data, although each quantity apparently follows a slightly non-gaussian distribution (Figure 11).

This figure provides several interesting observations. It demonstrates a clear correlation structure between $\Delta\sigma_{nn}$, $\Delta T_G$, and $\Delta T_D$, which are all positively correlated (with correlation values generally larger than 0.5, see also Figure S5 in Supporting Information S1). The most correlated quantities are $\Delta T_G$ and $\Delta T_D$, as inferred in the previous section, but it is also clear that both gouge and damage thickness variations are correlated with normal stress variations on the fault. Based on these joint distributions, we can derive some first order linear relationships between these three quantities (derived from their correlation matrix), which express the general trends existing between local deviations to the spatial means. We observe, for example, that, for $\sigma_3 = 90$ MPa, a local deviation in damage thickness $\Delta T_D$ is likely to be close to 2.16 times the local deviation in gouge thickness $\Delta T_G$.

Figure 13 allows to quantify the ratios associated with the first order linear relationships between $\Delta\sigma_{nn}$, $\Delta T_G$, and $\Delta T_D$, for the three confining stresses. As $\sigma_3$ increases, the ratio $\Delta T_D/\Delta\sigma_{nn}$ increases as well, which indicates that a certain local deviation in normal stress statistically induces a larger deviation in the damage thickness. This effect cumulates with the fact that the deviations in normal stress are much larger as $\sigma_3$ increases, and leads to the observation of Figure 10b, where the standard deviation of $T_D$ is highly sensitive to the confining stress. Conversely, the ratio $\Delta T_G/\Delta\sigma_{nn}$ essentially decreases with $\sigma_3$, which leads to a large increase of the ratio $\Delta T_D/\Delta T_G$. Hence, under a low level of confinement ($\sigma_3 = 45$ MPa), a given deviation to the average damage thickness correlates with a roughly equal deviation in gouge thickness ($\Delta T_D/\Delta T_G = 0.97$); while at large confining stresses ($\sigma_3 = 180$ MPa), the relative deviation in gouge thickness is about three times smaller ($\Delta T_D/\Delta T_G = 3.13$). This confirms the attenuating role of the gouge layer regarding these spatial heterogeneities.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** First order linear relationships between $\Delta\sigma_{nn}$, $\Delta T_G$, and $\Delta T_D$ (note that $\Delta\sigma_{nn}$ is expressed in MPa, and $\Delta T_G$, and $\Delta T_D$ are expressed in $\mu$m).</figcaption></figure>

Some additional information regarding the spatio-temporal structure of these correlations is brought by Figure 14 for the case M90. The blue curves correspond to the average autocovariograms of $\sigma_{nn}$, $T_G$, and $T_D$ during the last millisecond of the simulation, identified as a period of statistical steady state (in the sense that the statistical properties and spatial correlation structures of the observed heterogeneities stabilize). These curves show a typical decrease of the autocorrelation function as the lag distance is increased: Starting from a perfect correlation of 1 at a lag distance of 0, it decreases to much lower values after a certain distance. The lag distance at which the autocorrelation drops below $e^{-1}$ is conventionally termed “autocorrelation length”, and characterizes the size of the typical structures observed in a random signal. In the case of $\sigma_{nn}$, the autocorrelation length is found equal to 650 $\mu$m, while it is equal to 290 and 230 $\mu$m for $T_G$, and $T_D$, respectively. This indicates that the typical regions of high or low values are more extended for the stress profiles than for the damage and gouge thickness profiles. A second point of interest of these autocovariograms is the existence of a secondary peak of positive correlation, which indicates a certain level of spatial periodicity. It is interesting to notice that, in the case M90, these peaks occur at a similar distance for the three quantities: 3.30, 3.70, and 3.32 $\mu$m for $\sigma_{nn}$, $T_G$, and $T_D$, respectively. This can be interpreted as an identical periodicity in the three signals.

To get a sense of temporal evolution of this correlation structures, Figure 14 also provides maps of the evolution of these autocorrelation functions in time. This is completed by the curves providing the temporal evolution of the correlation coefficient between $\Delta\sigma_{nn}$, $\Delta T_G$, and $\Delta T_D$ (taken pairwise), in the left-hand part of the figure. At the beginning of the simulation (i.e., during the elastic loading and until 1 ms after the first major stress drop), we observe that the autocorrelation function of $\sigma_{nn}$ has a very sharp initial decrease (with a low autocorrelation length of ~200 $\mu$m), and no clear secondary peak: in that state, the profiles of normal stress are therefore disordered, with no clear periodicity, and sharp fluctuations in space. This observation also holds for $\Delta T_G$, and $\Delta T_D$, and is confirmed by the fact that the levels of spatial correlation between the normal stress and the gouge/damage

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Spatial correlation structure of $\sigma_{nn}$, $T_G$ and $T_D$, for a confining stress of 90 MPa. From left to right: (i) Evolution in time of the correlation coefficient of each pair of quantities; (ii) Evolution in time of the autocovariogram of $\sigma_{nn}$ and curve of the average autocorrelation function of $\sigma_{nn}$ during the last millisecond of simulation (in blue); (iii) Same plots for $T_G$; (iv) Same plots for $T_D$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Spatial correlation structure for different confining stresses and for the three quantities of interest: (a) Autocorrelation length; (b) Distance to the first major autocorrelation peak. Note that all data points correspond to averages on the last millisecond of simulation. Minor autocorrelation peaks are marked with hollow circles.</figcaption></figure>

thicknesses remains low ($\rho_{\sigma_{nn}-T_G}$ and $\rho_{\sigma_{nn}-T_D}$ are typically close to 0.2). The final, stable correlation structure is therefore not established on the fault at that stage. However, as seismic cycles keep developing and as the sliding distance increases, the profiles of $\sigma_{nn}$, $T_G$, and $T_D$ get more and more structured: between 1.5 and 4 ms, the correlation levels between stress and fault damage increase from 0.2 to 0.6, the autocorrelation distances of $\sigma_{nn}$, $T_G$, and $T_D$ increase toward their final values, and the secondary peaks progressively develop. After $t = 4$ ms, these quantities do not seem to evolve any more, and the steady state pictured in the blue curves is reached. It should be noted that, based on the maps of Figure 14, it is not possible to derive a clear causality relationship between the different quantities: variations in $\sigma_{nn}$ may be the cause of variations in $T_G$, and $T_D$, or it could be the opposite. What is clear, however, is that these variations evolve together and get more and more synchronized in space as the fault matures.

Cases M45 and M90 generally exhibit similar results (progressive increase of the correlations $\rho_{\sigma_{nn}-T_G}$ and $\rho_{\sigma_{nn}-T_D}$, of the autocorrelation lengths of $\sigma_{nn}$, $T_G$, and $T_D$, and appearance of secondary autocorrelation peaks, Figure S6 in Supporting Information S1), although trends are a bit less clear. These results are summarized in Figure 15 as a function of $\sigma_3$, considering the last millisecond of each simulation. Figure 15a shows that the autocorrelation lengths of $T_G$ and $T_D$ are always close to each other (typically between 200 and 300 µm), irrespective of the value of $\sigma_3$. The autocorrelation length of $\sigma_{nn}$ is systematically larger (i.e., between 400 and 650 µm), which indicates less sharp variations in space. Figure 15b provides the position of the secondary peaks of autocorrelation, and shows that the spatial synchronization between $\sigma_{nn}$, $T_G$, and $T_D$ is valid in all cases. The spatial period seems to be larger in the intermediate M90 case, but the limited number of simulations does not allow to be conclusive on the general relationship between this period and $\sigma_3$. Minor correlation peaks extracted from Figure 14 and Figure S6 in Supporting Information S1 are also provided, and seem to cluster around a period of ~1 mm. This corresponds to the typical spatial periods observed in the effective friction profiles (Figure 6 and Figure S3 in Supporting Information S1), and suggests an underlying spatial structure of the fault state (damage and gouge thickness) which would be associated with friction rather than normal stress.

To complement this analysis, we analyze the Power Spectral Densities (PSD) of the unfiltered stress profiles. Four such PSD are plotted in Figures 16a–16d, for both the normal ($\sigma_{nn}$) and the tangential ($\sigma_{nt}$) stress profiles on the footwall and on the hanging wall, in the final state of simulation M90. These PSD all exhibit a first drop followed by a clear peak close to 300 m⁻¹, which corresponds to the large-scale periodicity emphasized in the previous paragraphs. Beyond this peak, a power-law decay is observed. Beyond a wavenumber of ~6,000 m⁻¹, this decay becomes blurred by the noise associated with the grain-scale contact algorithm at the boundary between the degradable and continuum regions of the model. We therefore extract a decay exponent $-\alpha$ (with $\alpha > 0$) in the
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Spectral analysis of the stress profiles. (a) Power Spectral Density (PSD) of the profile of normal stress $\sigma_{nn}$ on the footwall at the end of simulation M90 (a power law fit is performed between 500 and 4,000 m$^{-1}$, shown as a black segment in the gray region); (b) PSD of $\sigma_{nn}$ on the hanging wall; (c) PSD of the profile of tangential stress $\sigma_{nt}$ on the footwall; (d) PSD of $\sigma_{nt}$ on the hanging wall; (e) Evolution of the four PSD exponents $\alpha$ of $\sigma_{nn}$ and $\sigma_{nt}$ on both walls along time for the simulation M90 (colors identical to (a–d.); (f) Averaged and smoothed PSD exponent $\alpha$ along time for the simulations M45, M90, and M180.</figcaption></figure>

range 500–4,000 m$^{-1}$, corresponding to spatial periods of 250 µm (close to the thickness of the damaged layer) to 2 mm (close to the large-scale periodicity of the fault heterogeneity), marked in gray in Figure 16. The evolution of $\alpha$ in time is plotted in Figure 11e for the four evolving stress profiles. Since the range of wavenumbers used for the extraction of $\alpha$ is limited, this estimate is subjected to a high level of noise, but the curves of Figure 11e nevertheless exhibit a clear and systematic evolution. The PSD exponent of normal stress is initially close to a value of 2.6, while that of shear stress is initially much lower, and close to 1.4. This observation holds for other values of the confining stress, see Figure S7 in Supporting Information S1. However, during the elastic loading, the four PSD exponents seemingly converge until they reach a similar value close to 2.0, and then follow a common evolution. This convergence is particularly clear in the M180 case (Figure S7 in Supporting Information S1). $\alpha$ then steadily increases up to values close to 3.0 as seismic cycles accumulate, and then decreases until reaching a plateau close to 2.2, until the end of the simulation. In Figure 16f, the four curves are averaged in one single curve and smoothed in order to ease the reading, and this averaged $\alpha$ is plotted for the three confining pressures. The general evolution previously described is also valid for the cases M45 and M180, although with some variations and with a clear time delay (associated with the duration of the initial elastic loading). We checked that the fitting of a power law on these PSD at any point in time remained legitimate on the same range of wavenumbers, and it appears to be the case. A point of caution, however, is that the locus of stress measurements (i.e., the boundary of the continuum region) is fixed in space while the damage front is progressing toward it. While this could be the cause of an artificial roughening of the measured stress field, we do not observe any major difference between the cases M45 and M180 (Figure S1 in Supporting Information S1), while the damage front gets much closer in the latter case. It therefore seems likely that the trends depicted in Figure 16f remain meaningful.
# 6. Discussion

## 6.1. Persistent Stress Heterogeneities

The numerical results reported in the previous sections provide new insights into the emergence and characteristics of fault complexity. Despite a fault zone with initially homogeneous mechanical properties, heterogeneous stress patterns start to appear early in the simulation, even before the first major stress drop (Figure 2). They then develop in time as the sliding distance increases. The stress variability is much larger in space (variations in fault local properties and structure) than in time (progressive loading and sudden unloading associated with seismic cycles), as quantified in Figure 7. It indicates that the stress patterns are essentially stationary on the fault, despite the numerous sliding events occurring at each location. This is in line with various observations reported in the field and in the lab at various scales (Ben-Zion & Dresen, 2022). To cite a few examples, Aubry et al. (2018) reported stress concentrations (with a size of the order of ~100 µm), apparently active during several stick-slip cycles, evidenced by localized shear heating in a lab-scale fault; Cebry et al. (2022) identified regions of repeating small events on a meter-scale lab fault, and identified them as persistent secondary asperities (with a size of several cm), developing from small heterogeneities in gouge distribution; Yamashita et al. (2021) provided direct measurements of heterogeneous shear stress profiles along a meter-scale lab fault, and confirmed that the stress distribution remains essentially unaffected by the large number of sliding events taking place in the interface; At the regional scale, Schmittbuhl et al. (2006) were able to infer the normal stress field on the Nojima fault before and after the 1995 Kobe earthquake, and found that this field was both highly variable in space and not much affected by the earthquake, which allowed them to conclude that the stress distribution and the asperities (about 10 km in size) were intrinsic properties of the fault itself; This comes in addition to the numerical study reported in Radiguet et al. (2013), which showed that shear stress concentrations “survive” the passing of several rupture fronts, although it was reported in the same study that a certain stress concentration only kept 1.6% of its initial value after the propagation of five slip events (presumably because no permanent fault feature, either structural or geometrical, was associated with this stress concentration).

## 6.2. Stress Probabilistic Distribution and Resulting Effective Friction

Profiles of shear stress and of normal stress were found to be strongly positively correlated in space and time (Figure S2 in Supporting Information S1), and the heterogeneities gained in amplitude as the sliding distance increased (Figure 7). Plotting the joint distribution of normal and shear stresses as a 2D histogram, we observe that the distributions strongly differ from the expected average loading path that would exist on the interface if it was perfectly homogeneous. This means that this heterogeneity is not related to some lack of synchronization of a similar seismic cycle between different regions of the fault, but rather to the local structure of the fault, which locally modifies this loading path (Figure 5). A possible explanation would be related to the rotation of the stress tensor either in the damaged region (Faulkner et al., 2006) or in the gouge layer (Giorgetti et al., 2019), although this seems unlikely since such a rotation is not expected to modify the normal and shear stresses resolved on the fault plane. Normal stress was found to follow a non-normal probability density function, with a somewhat triangular, bonded, essentially symmetric shape (Figure 11). Effective friction was found to generally vary between 0.1 and 0.3 for all confining stresses, and friction patterns in space were observed to be just as persistent and as structured as the stress patterns, albeit with an apparently different periodicity (Figure 6 and Figure S3 in Supporting Information S1).

## 6.3. Heterogeneous Damaging Process

An interesting feature of the present numerical model is the possibility for the fault rock to damage and to produce a layer of granular gouge. We indeed observe in all simulations that a damage profile develops, with a damage gradient evolving from the gouge layer (full damage) to the intact rock (no damage), as commonly observed in the field and in the lab (see for example Faulkner et al., 2006; Mitchell & Faulkner, 2009; Noël et al., 2023). We show that patterns in gouge and damage thickness also appear along the fault, and develop during sliding. Just like stress patterns, their locations are persistent in time (Figure 9 and Figure S4 in Supporting Information S1), with the major difference that damage is a cumulative process since no healing is implemented in the simulations: these patterns therefore come in addition to the continuous growth of the gouge layer and of the damaged region. Both gouge and damage thicknesses increase with slip, but at a decreasing rate, which is in good agreement with both lab (Boneh et al., 2014) and field (Sagy & Brodsky, 2009) estimates.
In Kwiatek et al. (2024), Westerley granite samples with a rough fault (spontaneously nucleated from notches during the test) were submitted to loading conditions similar to that used in the present model. Among other findings, many local slip events attributed to asperity failures were measured, and the evolution of the aftershock rates and of the b-values associated with the resulting catalog suggested a general smoothing of the fault surface. The authors postulated that fault damage and the evolution of roughness and stress across multiple stick-slip cycles tends to smooth small-scale asperities but generates persisting large-scale topography. It should be noted however that these results were obtained in the case of a fault which is initially much rougher than ours. As reported in Wang et al. (2024), starting with a planar and smooth fault also leads to the build-up of a thin layer of gouge, and triggers a certain roughening (see also the results from Aubry et al. (2020) on marble). In Wang et al. (2024), stress heterogeneity and roughness were found to control hypocenter distribution, and events clustering was shown to be closely correlated with the slip-parallel gradient of surface topography, with estimated source sizes spanning from 0.1 to 2.6 mm for acoustic events, not unlike the present findings. This is in overall agreement with the suggestion by Ben-Zion and Sammis (2003), that after a sufficient number of seismic cycles faults locally reach a critical state of disorder. An interesting consequence of this progressive maturation of the fault plane toward an established state was reported in Perrin et al. (2016), where the regions of largest slip of large earthquakes were shown to be highly correlated with the level and gradient in fault maturity. In future studies, such a prediction could be tested in our model by analyzing the slip profiles of large events and comparing them to the heterogeneous levels of damage and gouge thickness we reported in the present work.

## 6.4. Correlations Between Stress, Damaging, and Gouge Production

The variations in normal stress, in gouge thickness, and in damage thickness (expressed as a deviation to the instantaneous spatial mean to correct from temporal trends) were found to be strongly positively correlated (Figure 12 and Figure S5 in Supporting Information S1). In other words, regions of larger normal and shear stress (which are generally persistent in time) also are regions of increased gouge thickness and of increased damage intensity. Gouge accumulation is observed to be correlated with damage, but with a certain amount of smoothing in space, which is not surprising: gouge production is the direct consequence of rock damage, but the larger mobility of gouge grains allows them to spread along the fault and to attenuate the evolution and the acceleration of damage. Such a correlation between gouge thickness and fault geometry was also inferred in the field in Sagy & Brodsky, 2009, on a normal fault with a cumulated slip of 100–200 m in an andesite/sandstone setting. In this study, a feedback cycle was proposed: (a) slip on the fault produces a layer of wear particles, which has regions of variable thickness; (b) these variations of thickness induce internal deformation of the granular non-cohesive layer; (c) this leads to a contrast in rheology which localizes slip, and therefore injects more new particles into the granular layer; (d) eventually, this process generates fault topography. Our model cannot reproduce the full complexity of a real fault zone: it has accumulated a much lower amount of slip, and the complex granular processes in the active layer (e.g., comminution, healing, etc.) are not included in the simulations. However, the conclusion of Sagy & Brodsky, 2009 stating that asperities are both geometrical and rheological inhomogeneities is supported by the present results, and our model informs on the very early stages of such a long-term process.

Also interesting is the correlation between damage and normal stress. This could be interpreted in two different ways: a local increase in normal stress could lead to more damage in the same area; or a locally higher damage could lead to larger normal stress through damage-related restrained dilation. Both causalities could also be valid at the same time. In any case, this observation is reminiscent of a field study which showed that the intensity of shear heating (evidenced by the measured thickness of pseudotachylyte and possibly related to a locally larger normal stress) was correlated with a larger fragmentation of the wall rock (Petley-Ragan et al., 2019).

## 6.5. Fault Evolution

Before the deviatoric loading, the fault is mechanically homogeneous and planar, and the only heterogeneity is related to the random shapes of the grains. However, as illustrated in Figures 2, 6, and 9, Figures S1, S3, and S4 in Supporting Information S1, stress heterogeneities start to appear well-before the first major stress drop (either because of the random microstructure of the discrete region or because of the initial irregularity of the surfaces), in conjunction with localized (not pervasive) damaging and gouge production related to minor sliding events (described in Mollon et al., 2023). Main structures however only appear after this major event, which establishes a continuous layer of gouge in the interface and initiates the seismic cycles to follow. From thereon, the heterogeneity characteristics of different quantities follow different histories. The variability in normal stress is
generally found to mildly increase with slip, although trends are not clearly identical for different confining stresses (Figure 7). The variability in damage thickness increases steadily during the whole simulation (Figure 10b), without modifying the damaging spatial patterns (Figure 9 and Figure S4 in Supporting Information S1): areas with a large damage are prone to damage faster upon further sliding. Provided that damage is positively correlated with local normal stress, this observation points toward a positive feedback loop between local normal stress and damage rate (considered as deviations to the spatial mean). In contrast, both fault roughness and variability in gouge thickness seem to reach a rather constant level after the first major stress drop (despite the progressive growth of the gouge layer). Hence, regarding the gouge layer and its geometrical boundaries, the fault very quickly reaches a state of established disorder (Figures 10a and 10c), despite the increasing heterogeneities in stress and damage. In Noël et al. (2023), laboratory sliding tests were performed on bare rock surfaces of highly porous sandstone (among others) under confining stresses of 20 MPa and up to a final sliding distance of ~25 mm. Sandstone was intensely abraded during sliding: the interface progressively got filled with a layer of quartz grains pulled out from the surface and then submitted to comminution, until the end of the tests where the bare surfaces were separated by a layer of submicrometric gouge grains about 100–200 µm-thick. Rock walls exhibited a damage level that decreased with the distance to the fault plane. These observations are qualitatively consistent with the predictions of our model. An interesting point is that these experiments demonstrated that, as slip accumulates in this gouge layer, stability properties of the fault might evolve. A decrease of the Rate-and-State Friction parameters (a–b) and $D_c$ (both rendering the fault more prone to unstable sliding) was systematically observed, in particular during the initial rapid build up of the gouge layer. This was however related to the development of localization patterns (R1-bands and Y-bands, see Casas et al., 2023), and these are features that cannot be captured in our model because our grains are too large and we do not simulate comminution. Our findings however enrich this picture by indicating that the rate of gouge thickening is variable in space, and that the local evolution of the fault toward instability depicted in Noël et al. (2023) is certainly prone to the same spatial variability.

Figure 14 and Figure S6 in Supporting Information S1 show that the correlation between gouge and damage thickness is well-established from the first major stress drop. In contrast, the correlation between damage and normal stress steadily increases from ~0.2–0.5 during the first seismic cycles to ~0.6–0.7 after a sliding of ~400 µm. This result agrees with the observations from Cebry et al. (2022), where experiments of laboratory earthquake in a quartz gouge layer embedded in an elastic medium are reported. Despite the initial homogeneity of the gouge layer and of the elastic medium, these authors describe the appearance of so-called “secondary asperities”, which are reported to persist over many stick-slip cycles. They initiate from small heterogeneities in the initial gouge distribution and affect the local slip dynamics (peak velocity, seismic radiations, etc.), with more and more influence as cumulated sliding increases. This indicates that these asperities develop naturally as part of the evolving fault fabric and stress redistribution, and grow stronger throughout experiment. This observation agrees well with our predictions of persistent asperities materialized by increased stress. Besides, since the experiments reported in Cebry et al. (2022) do not imply any fault rock damage nor gouge production, we may infer that the primary cause for the appearance of these asperities stem from purely granular, fabric-related, phenomena. Our model brings the additional conclusion that, once initiated, these heterogeneities grow though a positive feedback though fault rock damage, additional gouge production, and local increase of the stress level.

## 6.6. Spatial Structure

The analysis of the spatial autocorrelation functions of normal stress, damage thickness and gouge thickness indicate that these quantities obey a clear spatial structuration (Figure 14 and Figure S6 in Supporting Information S1). This structuration, however does not appear immediately during the first stress drop but rather needs a large number of seismic cycles, corresponding to a slip of ~400 µm, to completely establish. This amount of slip seems to be also required to reach a maximum correlation between damage and stress (see previous paragraph). From thereon, a statistical steady state seems to be reached, where stress and damage spatial distributions have stabilized their spatial structure and tuned to maximize their correlation. In this established state, we observe that the autocorrelation length (i.e., the typical size of notable features in a random field) is larger for the normal stress (400–650 µm depending on the confining stress) than for the damage and gouge thicknesses (200–300 µm). However, the secondary autocorrelation peak indicates that the three quantities are submitted to a large amount of periodicity, with very similar periods in a given simulation (ranging between 2 and 3.5 mm depending on $σ_3$). Combining these two observations with the strong level of correlation between these quantities, we can confirm
that damage and stress are strongly synchronized in space, but that the damage distribution is submitted to sharper
spatial variations. The fact that these variations are similar but somewhat smoothed for the stress field could
indicate that stress is the consequence while damage is the cause, but data do not allow to be conclusive on this
point. Faults tend to an established quadratic roughness close to 8 µm for all confining stresses, while the
periodicity of the gouge thickness profile (closely related to the fault roughness) is close to 3 mm. This leads to an
estimate of a typical height-to-wavelength ratio of ~0.003 for the typical geometric asperity in the statistical
steady state of the fault. This estimate is in good agreement with field data provided in Sagy & Brodsky, 2009
(0.005) or in Power and Tullis (1991) (0.001–0.01). It is also interesting to mention that the asperities reported
through thermal measurements in the experiments of Aubry et al. (2018) (from which the present model is
inspired) had a typical size of 100 µm and a typical period of 450 µm. Although these asperities were measured
and characterized in the direction orthogonal to sliding, these measurements are not inconsistent with the present
numerical results.

At a scale below that of the main periodicity described in the previous paragraph, spatial structure is still present in the stress field under the form of a power-law decay of the PSD of stress fields, characterized by a time-varying α exponent. This hints at a possible self-affinity of the stress fields that spontaneously develop in our model, although the range of wavenumbers on which a power-law can be fitted is not large enough to be conclusive.

Previous field (Candela et al., 2012), lab (Aghababaei et al., 2022; Persson, 2014) and numerical (Milanese et al., 2019) works in fault mechanics and tribology demonstrated that the topography of rubbed surfaces tends to a steady state of complexity, with a constant roughness and a constant Hurst exponent. This constant value, which is close to 0.6 in the slip-parallel direction of faults, is believed to be the product of a dynamic equilibrium between roughening and smoothening processes. The self-affine character of the stress field of an idealized fault was demonstrated in Pétrélis et al., 2023 to control the b-value of its earthquake cycles, and such a correspondence will be investigated in future work on our numerical database.

6.7. Influence of the Confining Stress

Numerical results show that the confining stress has a limited effect on a number of phenomena related to fault heterogeneity. The gouge thickness and its spatial variability, for example, are only moderately influenced by σ₃. This is also the case of the fault roughness, and more generally of the spatial structuration of the fault (including periodicities, autocorrelation lengths and correlations among variables, Figure 15). Some parameters, however, are more strongly confinement-dependant. This is the case of the variability in stress and friction, which seems to be relatively smaller as σ₃ increases (Figures 5 and 11). This implies a certain (limited) homogenization of the fault under larger mean stresses. In contrast, both the average damage rate and the amplitude of the damage thickness heterogeneity increase monotonically with σ₃ (Figures 10b and 11). Hence, when the confining stress increases, a given deviation in normal stress leads to a larger deviation in damage thickness. Gouge thickness appears to be much less sensitive thanks to granular reorganization. This draws the picture of faults which get more mechanically homogeneous (normal stresses), more structurally heterogeneous (damage thickness), and geometrically identical (roughness and gouge thickness, spatial structuration) as the confining stress increases.

7. Conclusion

The numerical results presented in this paper provide a new perspective on the emergence, evolution, and
characteristics of fault heterogeneity, and reveal strong positive correlations between mechanical (stress fields)
and structural (host rock damage intensity and fault gouge thickness) fault properties. Gouge production distri-
bution is a direct consequence of rock damage distribution (with a spatial variability attenuated by grain-scale
mobility), which in turn appears to develop a mutual positive feedback with the normal stress profile on the
fault: regions of larger normal stress are submitted to a higher damage rate, which in turn promotes a local in-
crease of the normal stress (e.g., through damage-related restrained dilation). A clear spatial structure emerges
after a certain amount of slip on the fault (and a certain number of sliding event), with typical feature sizes and
spatial periodicities seemingly independent on the confining stress. A key question is the origin of the spatial
scales that spontaneously appear in this simplified system. If they are essentially driven by granular phenomena
occurring within the gouge layer, they may be related to the grain size that was arbitrarily chosen in the model; but
if they are related to elasto-frictional length-scales associated with the dynamics of sliding events, they could
prove to be intrinsic fault characteristics. This will be the topic of future studies.
An unexpected point is that the spatial structure of the effective friction (Figure 6 and Figure S3 in Supporting Information S1) does not seem to follow that of the geometric asperities of the fault, but rather to obey to a different (shorter) spatial periodicity for all values of $\sigma_3$. The causes for this discrepancy remain unclear based on the present data. What is clear, however, is that local effective friction is essentially uncorrelated with the local damage level, gouge thickness, and normal stress (Figure 5). In future works, we will keep investigating this data set by focusing on the sliding kinematics during the numerous sliding events occurring in each simulation. An important work to be done is an accurate cataloging of all events in a given simulation (epicenter, size, duration, slip, velocity, stress drop, etc.). An important line of investigation might be to determine what is the feedback between fault heterogeneity and slip and stress drop distributions, as well as rupture nucleation, initiation, propagation and arrest. The consequences of the structuration of the fault on the seismic cycles, and their evolution during the maturation of the fault toward its established state of disorder, will be at the core of our future research.

## Data Availability Statement

All simulations were performed with the open-source software MELODY version 3.94, developed by the first author and described in (Mollon, 2018). The data generated is available at (Mollon et al., 2023).

## Acknowledgments

We thank the Associate Editor, as well as G. Dresen and two other anonymous reviewers for allowing us to improve our manuscript. 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

Aghababaei, R., Brodsky, E. E., Molinari, J.-F., & Srinivasan, C. (2022). How roughness emerges on natural and engineered surfaces. *MRS Bulletin*, 47(12), 1229–1236. https://doi.org/10.1557/s43577-022-00469-1

Agosta, F., & Aydin, A. (2006). Architecture and deformation mechanism of a basin-bounding normal fault in Mesozoic platform carbonates, central Italy. *Journal of Structural Geology*, 28(8), 1445–1467. https://doi.org/10.1016/j.jsg.2006.04.006

Albertini, G., Simon, K., Grigoriu, M. D., & Kammer, D. S. (2021). Stochastic properties of static friction. *Journal of the Mechanics and Physics of Solids*, 147, 104242. https://doi.org/10.1016/j.jmps.2020.104242

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

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(12), 274–282. 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

Bedford, J. D., Faulkner, D. R., & Lapusta, N. (2022). Fault rock heterogeneity can produce fault weakness and reduce fault stability. *Nature Communications*, 13(1), 326. https://doi.org/10.1038/s41467-022-27998-2

Ben-Zion, Y., & Dresen, G. (2022). A synthesis of fracture, friction and damage processes in earthquake rupture zones. *Pure and Applied Geophysics*, 179(12), 4323–4339. https://doi.org/10.1007/s00024-022-03168-9

Ben-Zion, Y., & Sammis, C. G. (2003). Characterization of fault zones. *Pure and Applied Geophysics*, 160(3–4), 677–715. https://doi.org/10.1007/PL00012554

Beroza, G. C., & Mikumo, T. (1996). Short slip duration in dynamic rupture in the presence of heterogeneous fault properties. *Journal of Geophysical Research*, 101(B10), 22449–22460. https://doi.org/10.1029/96JB02291

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

Bruhat, L., Klinger, Y., Vallage, A., & Dunham, E. M. (2020). Influence of fault roughness on surface displacement: From numerical simulations to coseismic slip distributions. *Geophysical Journal International*, 220(3), 1857–1877. https://doi.org/10.1093/gji/ggz545

Candela, T., Renard, F., Klinger, Y., Mair, K., Schmittbuhl, J., & Brodsky, E. E. (2012). Roughness of fault surfaces over nine decades of length scales: fault roughness. *Journal of Geophysical Research*, 117(B8), B08409. https://doi.org/10.1029/2011JB009041

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

Casas, N., Mollon, G., & Daouadji, A. (2023). Influence of grain-scale properties on localization patterns and slip weakening within dense granular fault gouges. *Journal of Geophysical Research: Solid Earth*, 128(3), e2022JB025666. https://doi.org/10.1029/2022JB025666

Cattania, C., & Segall, P. (2021). Precursory slow slip and foreshocks on rough faults. *Journal of Geophysical Research: Solid Earth*, 126(4), e2020JB020430. https://doi.org/10.1029/2020JB020430

Cebry, S. B. L., Ke, C. Y., Shreedharan, S., Marone, C., Kammer, D. S., & McLaskey, G. C. (2022). Creep fronts and complexity in laboratory earthquake sequences illuminate delayed earthquake triggering. *Nature Communications*, 13(1), 1–9. https://doi.org/10.1038/s41467-022-34397-0

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

Chester, F. M., Evans, J. P., & Biegel, R. L. (1993). Internal structure and weakening mechanisms of the san andreas fault. *Journal of Geophysical Research*, 98(B1), 771–786. https://doi.org/10.1029/92JB01866

Collettini, C., Niemeijer, A., Viti, C., & Marone, C. (2009). Fault zone fabric and fault weakness. *Nature*, 462(7275), 907–910. https://doi.org/10.1038/nature08585
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

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

Dresen, G., Kwiatek, G., Goebel, T., & Ben-Zion, Y. (2020). Seismic and aseismic preparatory processes before large stick-slip failure. *Pure and Applied Geophysics*, **177**(12), 5741–5760. https://doi.org/10.1007/s00024-020-02605-x

Dublanchet, P., Bernard, P., & Favreau, P. (2013). Interactions and triggering in a 3-D rate-and-state asperity model. *Journal of Geophysical Research: Solid Earth*, **118**(5), 2225–2245. https://doi.org/10.1002/jgrb.50187

Faulkner, D. R., Mitchell, T. M., Healy, D., & Heap, M. J. (2006). Slip on “weak” faults by the rotation of regional stress in the fracture damage zone. *Nature*, **444**(7121), 922–925. https://doi.org/10.1038/nature05353

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

Ghosh, A., Vidale, J. E., & Creager, K. C. (2012). Tremor asperities in the transition zone control evolution of slow earthquakes: Tremor asperities control slow quakes. *Journal of Geophysical Research*, **117**(B10), B10301. https://doi.org/10.1029/2012JB009249

Giorgetti, C., Tesei, T., Scuderi, M. M., & Collettini, C. (2019). Experimental insights into fault reactivation in gouge-filled fault zones. *Journal of Geophysical Research: Solid Earth*, **124**(4), 4189–4204. https://doi.org/10.1029/2018JB016813

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

Gounon, A., Latour, S., Letort, J., & El Arem, S. (2022). Rupture nucleation on a periodically heterogeneous interface. *Geophysical Research Letters*, **49**(20), e2021GL096816. https://doi.org/10.1029/2021GL096816

Guérin-Marthe, S., Kwiatek, G., Wang, L., Bonnelye, A., Martínez-Garzón, P., & Dresen, G. (2023). Preparatory slip in laboratory faults: Effects of roughness and load point velocity. *Journal of Geophysical Research: Solid Earth*, **128**(4), e2022JB025511. https://doi.org/10.1029/2022JB025511

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

Harbord, C. W. A., Nielsen, S. B., De Paola, N., & Holdsworth, R. E. (2017). Earthquake nucleation on rough faults. *Geology*, **45**(10), 931–934. https://doi.org/10.1130/G39181.1

Hu, J., Song, H., Sandfeld, S., Liu, X., & Wei, Y. (2022). Breakdown of Archard law due to transition of wear mechanism from plasticity to fracture. *Tribology International*, **173**, 107660. https://doi.org/10.1016/j.triboint.2022.107660

Kandula, N., Cordonnier, B., Boller, E., Weiss, J., Dysthe, D. K., & Renard, F. (2019). 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

Kemeny, J. M., & Hagaman, R. M. (1992). An asperity model to simulate rupture along heterogeneous fault surfaces. In C. G. Sammis, M. Saito, & G. C. P. King (Eds.), *Fractals and chaos in the Earth sciences. Pageoph topical volumes*. Birkhäuser. https://doi.org/10.1007/978-3-0348-6191-5_3

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

Kim, Y.-S., Peacock, D. C. P., & Sanderson, D. J. (2004). Fault damage zones. *Journal of Structural Geology*, **26**(3), 503–517. https://doi.org/10.1016/j.jsg.2003.08.002

Kwiatek, G., Martínez-Garzón, P., Goebel, T., Bohnhoff, M., Ben-Zion, Y., & Dresen, G. (2024). Intermittent criticality multi-scale processes leading to large slip events on rough laboratory faults. *JGR Solid Earth*, **129**(3), e2023JB028411. https://doi.org/10.1029/2023JB028411

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

Lebihain, M., Roch, T., Violay, M., & Molinari, J.-F. (2021). Earthquake nucleation along faults with heterogeneous weakening rate. *Geophysical Research Letters*, **48**(21), e2021GL094901. https://doi.org/10.1029/2021GL094901

Leeman, J. R., Saffer, D. M., Scuderi, M. M., & Marone, C. (2016). Laboratory observations of slow earthquakes and the spectrum of tectonic fault slip modes. *Nature Communications*, **7**(1), 11104. https://doi.org/10.1038/ncomms11104

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

Luo, Y., & Ampuero, J. P. (2018). Stability of faults with heterogeneous friction properties and effective normal stress. *Tectonophysics, Physics of Earthquake Rupture Propagation*, **733**, 257–272. https://doi.org/10.1016/j.tecto.2017.11.006

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), 2219. https://doi.org/10.1029/2001JB000516

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., Mohamed, A., & 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

Milanese, E., Brink, T., Aghababaei, R., & Jean-François Molinari, J. F. (2019). Emergence of self-affine surfaces during adhesive wear. *Nature Communications*, **10**(1), 1116. https://doi.org/10.1038/s41467-019-09127-8

Mitchell, T. M., & Faulkner, D. R. (2009). The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: A field study from the atacama fault system, northern Chile. *Journal of Structural Geology*, **31**(8), 802–816. https://doi.org/10.1016/j.jsg.2009.05.002

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 [Software]. *Computational Particle Mechanics*, **5**(4), 517–527. https://doi.org/10.5281/zenodo.4305614
Mollon, G., Aubry, J., & Schubnel, A. (2021). Simulating melting in seismic fault gouge. *Journal of Geophysical Research: Solid Earth*, 126(6).
https://doi.org/10.1029/2020jb021485

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

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

Noël, C., Giorgetti, C., Scuderì, M. M., Collettini, C., & Marone, C. (2023). The effect of shear displacement and wear on fault stability: Laboratory constraints. *Journal of Geophysical Research: Solid Earth*, 128(4), e2022JB026191. https://doi.org/10.1029/2022JB026191

Papachristos, E., Stefanou, I., & Sulem, J. (2023). A discrete elements study of the frictional behavior of fault gouges. *Journal of Geophysical Research: Solid Earth*, 128(1), e2022JB025209. https://doi.org/10.1029/2022JB025209

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

Perrin, C., Manighetti, I., Ampuero, J.-P., Cappa, F., & Gaudemer, Y. (2016). Location of largest earthquake slip and fast rupture controlled by along-strike change in fault structural maturity due to fault growth. *Journal of Geophysical Research: Solid Earth*, 121(5), 3666–3685. https://doi.org/10.1002/2015JB012671

Persson, B. N. J. (2014). On the fractal dimension of rough surfaces. *Tribology Letters*, 54(1), 99–106. https://doi.org/10.1007/s11249-014-0313-4

Petley-Ragan, A., Ben-Zion, Y., Austrheim, H., Ildefonse, B., Renard, F., & Jamtveit, B. (2019). Dynamic earthquake rupture in the lower crust. *Science Advances*, 5(7), eaaw0913. https://doi.org/10.1126/sciadv.aaw0913

Pétrélis, F., Chanard, K., Schubnel, A., & Hatano, T. (2023). Earthquake magnitude distribution and aftershocks: A statistical geometry explanation. *Physical Review E*, 107(3), 034132. https://doi.org/10.1103/PhysRevE.107.034132

Power, W. L., & Tullis, T. E. (1991). Euclidean and fractal models for the description of rock surface roughness. *Journal of Geophysical Research*, 96(B1), 415–424. https://doi.org/10.1029/90JB02107

Queener, C. A., Smith, T. C., & Mitchell, W. L. (1965). Transient wear of machine parts. *Wear*, 8(5), 391–400. https://doi.org/10.1016/0043-1648(65)90170-5

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

Reichelt, M., & Cappella, B. (2021). Large scale multi-parameter analysis of wear of self-mated 100Cr6 steel – a study of the validity of Archard’s law. *Tribology International*, 159, 106945. https://doi.org/10.1016/j.triboint.2021.106945

Renard, F., Weiss, J., Mathiesen, J., Ben-Zion, Y., Kandula, N., & Cordonnier, B. (1986). Critical evolution of damage toward system-size failure in crystalline rock. *JGR Solid Earth*, 123(2), 1969–1986. https://doi.org/10.1002/2017JB014964

Ripperger, J., Ampuero, J.-P., Mai, P. M., & Giardini, D. (2007). Earthquake source characteristics from dynamic rupture with constrained stochastic fault stress. *Journal of Geophysical Research*, 112(B4), B04311. https://doi.org/10.1029/2006JB004515

Romanet, P., Bhat, H. S., Jolivet, R., & Madariaga, R. (2018). Fast and slow slip events emerge due to fault geometrical complexity. *Geophysical Research Letters*, 45(10), 4809–4819. https://doi.org/10.1029/2018GL077579

Rubino, V., Lapusta, N., & Rosakis, A. J. (2022). Intermittent lab earthquakes in dynamically weakening fault gouge. *Nature*, 606(7916), 922–929. https://doi.org/10.1038/s41586-022-04749-3

Sagy, A., & Brodsky, E. E. (2009). Geometric and rheological asperities in an exposed fault zone. *Journal of Geophysical Research*, 114(B2), B02301. https://doi.org/10.1029/2008JB005701

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

Sammis, C. G., Nadeau, R. M., & Johnson, L. R. (1999). How strong is an asperity? *Journal of Geophysical Research*, 104(B5), 10609–10619. https://doi.org/10.1029/1999JB900006

Schmittbuhl, J., Chambon, G., Hansen, A., & Bouchon, M. (2006). Are stress distributions along faults the signature of asperity squeeze? *Geophysical Research Letters*, 33(13), L13307. https://doi.org/10.1029/2006GL025952

Scholz, C. H. (2019). *The mechanics of earthquakes and faulting* (3rd ed.). Cambridge University Press. https://doi.org/10.1017/9781316681473

Scott, D. R. (1996). Seismicity and stress rotation in a granular model of the brittle crust. *Nature*, 381(6583), 592–595. https://doi.org/10.1038/381592a0

Scuderì, M. M., Collettini, C., Viti, C., Tinti, E., & Marone, C. (2017). Evolution of shear fabric in granular fault gouge from stable sliding to stick slip and implications for fault slip mode. *Geology*, 45, 731–734. https://doi.org/10.1130/G39033.1

Selvadurai, P. A., & Glaser, S. D. (2017). Asperity generation and its relationship to seismicity on a planar fault: A laboratory simulation. *Geophysical Journal International*, 208(2), 1009–1025. https://doi.org/10.1093/gji/ggw439

Shipton, Z. K., & Cowie, P. A. (2001). Damage zone and slip-surface evolution over mm to km scales in high-porosity Navajo sandstone, Utah. *Journal of Structural Geology*, 23(12), 1825–1844. https://doi.org/10.1016/s0191-8141(01)00035-9

Sobolev, G. A., Ponomarev, A. V., Koltsov, A. V., & Smirnov, V. B. (1996). Simulation of triggered earthquakes in the laboratory, induced seismic events. *Pure and Applied Geophysics*, 147, 345–355.

Taboada, A., & Renouf, M. (2023). Rheology and breakdown energy of a shear zone undergoing flash heating in earthquake-like discrete element models. *Geophysical Journal International*, 233(2), 1492–1514. https://doi.org/10.1093/gji/ggad004

Tal, Y., Hager, B. H., & Ampuero, J. P. (2018). The effects of fault roughness on the earthquake nucleation process. *Journal of Geophysical Research: Solid Earth*, 123(1), 437–456. https://doi.org/10.1002/2017JB014746

Tesei, T., Collettini, C., Barchi, M. R., Carpenter, B. M., & Di Stefano, G. (2014). Heterogeneous strength and fault zone complexity of carbonate-bearing thrusts with possible implications for seismicity. *Earth and Planetary Science Letters*, 408, 307–318. https://doi.org/10.1016/j.epsl.2014.10.021

Wang, C., Elsworth, D., & Fang, Y. (2019). Ensemble shear strength, stability, and permeability of mixed mineralogy fault gouge recovered from 3D granular models. *Journal of Geophysical Research: Solid Earth*, 124(1), 425–441. https://doi.org/10.1029/2018JB016066

Wang, L., Kwiatek, G., Renard, F., Guérin-Marthe, S., Rybacki, E., Bohnhoff, M., et al. (2024). Fault roughness controls injection-induced seismicity. *Proceedings of the National Academy of Sciences of the United States of America*, 121(3), e2310039121. https://doi.org/10.1073/pnas.2310039121

Wibberley, C. A. J., Yielding, G., & Di Toro, G. (2008). Recent advances in the understanding of fault zone internal structure: A review. *Geological Society, London, Special Publications*, 299(1), 5–33. https://doi.org/10.1144/SP299.2

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
Yamashita, F., Fukuyama, E., Xu, S., Kawakata, H., Mizoguchi, K., & Takizawa, S. (2021). Two end-member earthquake preparations illuminated by foreshock activity on a meter-scale laboratory fault. *Nature Communications*, **12**(1), 4302. https://doi.org/10.1038/s41467-021-24625-4

Zielke, O., Galis, M., & Mai, P. M. (2017). Fault roughness and strength heterogeneity control earthquake size and stress drop. *Geophysical Research Letters*, **44**(2), 777–783. https://doi.org/10.1002/2016GL071700