# DEM Analyses of Cemented Granular Fault Gouges at the Onset of Seismic Sliding: Peak Strength, Development of Shear Zones and Kinematics

Nathalie Casas, Guilhem Mollon, Ali Daouadji

► To cite this version:

Nathalie Casas, Guilhem Mollon, Ali Daouadji. 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, 2022, 179 (2), pp.679-707. 10.1007/s00024-021-02934-5. hal-03659790

HAL Id: hal-03659790

https://hal.science/hal-03659790

Submitted on 12 May 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# DEM Analyses of Cemented Granular Fault Gouges at the Onset of Seismic Sliding: Peak Strength, Development of Shear Zones and Kinematics

NATHALIE CASAS,<fnref n="1" />,<fnref n="2" /> ID GUILHEM MOLLON,<fnref n="1" /> and ALI DAOUADJI<fnref n="2" />

**Abstract**—Fault zones usually present a granular gouge, com-ing from the wear material of previous slips. This layer contributes to friction stability and plays a key role in the way elastic energy is released during sliding. Considering a mature fault gouge with a varying amount of mineral cementation between particles, we aim to understand the influence of the strength of interparticle bonds on slip mechanisms by employing the discrete element method. We consider a direct shear model without fluid in 2D, based on a granular sample with angular and faceted grain shapes. Focusing on the physics of shear accommodation inside the granular gouge, we explore the effect of an increase of cementation on effective friction (i.e. stress ratio) within the fault. We find that brittleness and the overall shear strength are enhanced with cementation, especially for dense materials. For the investigated data range, three types of cemented material are highlighted: a poorly cemented material (Couette flow profile, no cohesion), a cemented material with aggregates of cemented particles changing the granular flow and acting on slip weakening mechanisms (Riedel shear bands R), and a highly-cemented material behaving as a brittle material (with several Riedel bands followed by fault-parallel shear-localization Y). Effective friction curves present double weakening shapes for dense samples with enough cementation. We find that the effective friction of a cemented fault cannot be directly predicted from Mohr-Coulomb criteria because of the heterogeneity of the stress state and kinematic constraints of the fault zone.

**Keywords:** Cemented gouges, faulting, shear zones, DEM, granular mechanics, tribology.

## 1. Introduction

During earthquakes, frictional sliding releases the stresses accumulated in the pre-stressed surrounding

medium. The fault gouge, identified as the wear material of previous slips, contributes to friction stability (Marone & Scholz, 1988) and plays an important role in the sudden energy release at the onset of seismic sliding (Sammis et al., 1987). The number and amplitudes of successive slips occurring within the gouge reduce the size of particles towards a fractal distribution, which also reduces pore spaces (Blenkinsop, 1991; Muto et al., 2015; Sammis & Biegel, 1989). This particle size distribution can be explained by grain breakage theory (Daouadji & Hicher, 2010). Within a mature fault gouge, mineral cementation coming from rock dissolution, melting, or other processes derived from previous slips, can fill remaining pore spaces between particles and change the global state of cohesion (Philit et al., 2018). (Lade & Overton, 1989) showed that, for low confining pressures, the increase of cementation and the associated tensile strength lead to an enhancement of effective friction (i.e. ratio of the tangential stress to the applied normal stress). This phenomenon gives birth to a new, stronger granular material combining its history, the state of initial density (i.e. porosity within the sample), and cementation (Schellart, 2000). During gouge shearing, at low confining pressure, the breakage of cement bridges between particles enhances local kinematic freedom and allows grains reorganization, such as the formation of shear plane failures and global dilation (Tengattini et al., 2014).

A wide range of gouge characteristics (mineralogy, thickness, particle size, friction...) and deformation conditions (pressure, slip velocity...) are believed to influence, if not control, a large part of slip mechanisms. Such gouge parameters have been studied in the literature from Lab or in-situ points of

<footnote n="1">¹ Univ. Lyon, INSA-Lyon, CNRS UMR5259, LaMcS, 69621 Villeurbanne, France. E-mail: nath27casas@gmail.com</footnote>

<footnote n="2">² Univ Lyon, INSA Lyon, GEOMAS, EA7495, 69621 Villeurbanne, France.</footnote>
view (Anthony & Marone, 2005; Biegel et al., 1989; Byerlee & Brace, 1968; Mair & Marone, 1999; Marone & Scholz, 1989; Sammis et al., 1987) and also numerically by the mean of discrete element modelling (DEM) (Cho et al., 2008; Da Cruz et al., 2005; Gao et al., 2018; Guo & Morgan, 2004; Morgan, 1999; Morgan & Boettcher, 1999; Zhao et al., 2012). Cementation and its influence on slip behaviour are less well understood in the field of fault mechanics. It may be due to the difficulty to introduce and control a cemented gouge material in experiments, but also to their atypical behaviour in introducing local force chains within the gouge. This is the reason why discrete element modelling (DEM) is thought to be a good way to complete knowledge on cemented materials and their influence on slip mechanisms within the gouge in terms of strength, slip weakening, and failure patterns. In DEM, the cementation can be schematized by cohesive and breakable bonds between particles. Bonded particle models are often used to represent cohesive laws within granular rocks (Cho et al., 2007; Potyondy & Cundall, 2004) or to represent cementation (Das et al., 2014; Estrada et al., 2010; Jiang et al., 2013; Tengattini et al., 2014). DEM studies have traditionally used circular or spherical shapes (Dorostkar et al., 2018; Guo & Morgan, 2004), for the sake of computational and conceptual simplicity. In the recent years however, a general trend is to introduce more and more complexity in these shapes, either to investigate the influence of specific shape features (e.g. elongation), or to improve the predictive character of DEM by making the synthetic grains shapes as close as possible to those of the targeted granular population, as used in the recent study of (Harmon et al., 2021). We place our study in this second framework, and try to reproduce the contour of particles observed in granular fault gouge with its local roughness, angularity, and faceting (Muto et al., 2015; Olgaard & Brace, 1983). The effect of angularity on shearing behaviour was already demonstrated in literature (Anthony & Marone, 2005), as well as the effect of surface roughness on shear localisation (Mollon et al., 2020).

Literature reports a large range of phenomena occurring during faulting, mainly observed under high pressure and velocity and within a very thin

slipping zone. The temperature rise can lead to melting (Di Toro et al., 2006; Niemeijer et al., 2011) or fluid pressurization (Rice, 2006). These heating processes are related to a large mechanical work rate within the slip zone (Di Toro et al., 2011) and are responsible for thermal weakening. The well-known critical slip distance $D_c$ is the slip distance corresponding to a decrease in frictional strength [and is part of rate and state friction laws (Marone, 1998)], and is in the order $10^{-5}-10^{-3}$ m for lab experiments at low slip velocity (Marone & Kilgore, 1993). Another contribution to weakening is related to fracture mechanics and slipping zone. The mechanism by which frictional resistance drops with sliding is known as “slip weakening” and may be the cause of unstable sliding, releasing energy within the system in the form of high slip rate and stress drop across the fault. It is thus influenced by both the roughness of the fault surface and the thickness of the gouge layer (Marone & Kilgore, 1993). The main slipping zone can contain cohesive or non-cohesive wear material (Rice & Cocco, 2002), and a micrometre-millimetre scale “slip distance” is observed. The granular literature related to fault mechanics often explores stick-slip instabilities, considering strain energy storage either at grain-scale (Dorostkar et al., 2018; Leeman et al., 2015) or in the loading system (Kasyap & Senetakis, 2021). Although such energy storage is not present in our simulation, it is to be kept in mind from the perspective of upscaling the frictional behaviours we report to actual fault systems. In this paper, we also restrict our focus on the most immediate (fracture-related) slip weakening phenomena and disregard all thermal effects, which are supposed to happen only after a certain amount of slip.

The objective of this paper is therefore to establish a link between the micro-mechanical and structural properties of a cemented gouge layer, and its rheological behaviour under shearing. Varying the percentage of cementation within the gouge, for different initial porosities, leads to a wide range of mechanical behaviours that can be compared to typically observed fracture patterns. Section 2 introduces the numerical method for direct shear experiment and the fault segment model in 2D ($2 \times 20 \text{ mm}^2$), involving two rough surfaces representing the rock walls separated by the cemented granular gouge. This
paper proposes to combine modern simulation tech-
niques for cemented granular material with DEM on
angular and faceted grains, as we can observe in real
granular fault gouges. Focusing on the strength of
contact bonds inside the granular medium, we
explore strength evolution, gouge kinematics, and
force chains within the gouge in a third section. The
effect of cementation and initial porosity on
mechanical behaviours and kinematics of shear bands
are also investigated. The last section offers a dis-
cussion on new insights and relations between
cementation within the gouge, shear localization, and
the relevance of Mohr–Coulomb theory for fault
models.

## 2. Numerical Framework and Sample Generation

### 2.1. Discrete Element Modelling with Angular and Faceted Grains

Since DEM was first proposed by (Cundall &
Strack, 1979), it has been applied several times to the
simulation of micro-scale behaviours inside fault
zones (Cho et al., 2008; Da Cruz et al., 2005; Gao
et al., 2018; Guo & Morgan, 2004; Morgan, 1999;
Morgan & Boettcher, 1999; Zhao et al., 2012). In this
numerical framework, each particle has its own
motion and trajectory, driven by Newton's laws of
motion. The motion is controlled by user-defined and
physics-based contact interactions and constitutive
laws. Since the contour of each grain is discretized by
a piecewise linear frontier with nodes and segments,
each contact considered in the code concerns a given
node from a grain A and a given segment from a
grain B (Fig. 1). From this node and this segment, we
can at any time compute a normal gap $\delta_n$ (obtained
by projecting the node on the segment) and a
tangential gap $\delta_t$ (integrated in time based on the
history of the relative motions of the node and the
segment in the tangential direction). $\delta_n$ can be either
negative (i.e. there is a small interpenetration
between the grains) or positive (i.e. there is a
separation distance between the grains). The contact
detection is bi-directional, meaning that each grain
can adopt the roles of grains A and B at the same
time. Edge-edge contacts hence do not need any

special treatment in this framework, since they are
automatically accounted for through to their extrem-
ity nodes. This contact algorithm is described in more
details in Mollon (2018a). To control overlapping, a
contact stiffness is introduced in the normal and
tangential directions ($k_n$ and $k_t$).

The code used for simulations, MELODY 2D
(Multibody Element-free Open code for Dynamic
simulation), is a C++ code allowing to simulate a
broad variety of granular media. In contrast with
conventional codes (circular particles, clump logic,
convex polygonal particles...), it can deal with any
arbitrary 2D shape, and behaviour of particles, from
rigid circular to highly compliant angular grains (this
latter case is outside of the strict DEM framework,
and uses a Multibody Meshfree Approach). New
results are expected for sheared granular gouges with
complex shapes of particles closer to observations
(angular and faceted grains). More details on soft-
ware are available in Mollon (2016, 2018b) and
details on DEM can be found in (Supporting Infor-
mation S1 and S2). The detailed equations used for
cemented material are presented in the next section.

## 2.2. Cementation as Bonded Particles Within the Gouge

Hereafter, a fault zone that has already been
sheared in previous slip episodes is considered, which
became mature enough to observe mineral cementa-
tion between particles. As the cement present a rather
“cohesive” behaviour (Riedmüller et al., 2001;
Wibberley et al., 2008), this cementation can be

<figure><img src="image_2.png" /><figcaption>Figure 1
Sketch of a typical contact between two grains A and B. $k_n$ and $k_t$ are respectively the normal stiffness and the tangential stiffness and $\delta_n$ and $\delta_t$ are the associated gaps used to control interpenetrations between particles</figcaption></figure>
explicitly simulated in DEM by considering breakable cohesive bonds between particles. A “Bonded Mohr–Coulomb” law adapted to the complex grain shapes is therefore applied. This contact law is close to the bonded-particle-model (BPM) from (Potyondy & Cundall, 2004) and presents two main statuses (intact or broken) described below:

(a) After compaction and before shearing, all contacts for which $\delta_n < (\delta_{\text{detection}} = C_{\text{num}}/k_n)$, where $k_n$ is the contact stiffness and $C_{\text{num}}$ is the numerical cohesion, receive the status “intact” (Fig. 2a).

If the contact is “intact”, the contact stresses are computed based on a purely elastic contact law:

$$ \sigma_n = k_n \delta_n \quad (1) $$

$$ \sigma_t = k_t \delta_t \quad (2) $$

If $[\sigma_n$ or $abs(\sigma_t)]$ exceeds the prescribed value of cohesion $C_{\text{num}}$, the status of the contact is updated to “broken” (Fig. 2b).

(b) If the contact is “broken”, (Fig. 2b), either because it is a former intact bond or because it is a new interparticle contact formed during particle reorganization, the following contact stresses are computed based on a purely frictional contact law:

$$ \text{if } (\delta_n > 0) \rightarrow \sigma_n = \sigma_t = 0 \quad (3) $$

$$ \text{if } (\delta_n < 0) \rightarrow \begin{cases} \sigma_n = k_n \delta_n & (4) \\ \sigma_t = \min(k_t \delta_t, \mu_{\text{num}} \sigma_n) & (5) \end{cases} $$

where $\mu_{\text{num}}$ is the inter-particle friction coefficient.

From these stresses, the associated contact forces (in the normal and tangential direction, as well as the associated torque) are computed on each grain, by considering that contact stresses act on a contact length $L_c$ (equal to the sum of half-lengths of the segments around contact nodes in grain A):

$$ F_n = L_c \sigma_n \quad (6) $$

$$ F_t = L_c \sigma_t \quad (7) $$

The numerical cohesion $C_{\text{num}}$ considered in the code cannot be directly related to real cementation in rocks and does not bear much physical signification by itself. A good way to associate a quantitative description to this numerical parameter is to define an associated level of cementation. In the initial state of our model, each cohesive bond between any pair of contacting grains requires a certain amount of mechanical energy for breaking. This energy is related to its tensile and tangential stiffnesses, its tensile and tangential strengths, and its contact length $L_c$. To clarify the total amount of energy that would be needed in order to break all the initial cohesive bonds, we normalize it with respect to a representative energy. This energy corresponds to a surface energy of 62 J/m² which was reported for the Chilhowee quartzite and considered as an upper limit for rock surface energy by Friedman et al., (1972). We apply this surface energy on the whole external surface of all grains present in the simulation (a unit length is considered in the third dimension for any

<figure><img src="image_3.png" /><figcaption>Figure 2</figcaption></figure>

a Intact bond. Initialization of the contact law. A cohesive law links all grains in contact. b Broken bond. When the force applied to the particles becomes higher than the cohesive strength, the bond is broken. The contact becomes cohesionless and follows a classical Mohr–Coulomb law with inter-particle friction only. Broken contacts cannot be cohesive again and this induces an augmentation of broken bonds during the shearing. Normal interpenetration is exaggerated in the figure for the sake of clarity
necessary purpose), and consider that this amount of energy corresponds to a complete cementation of the gouge layer. Hence, any initial state of cementation effectively introduced in the numerical samples by the means of the numerical cohesion $C_{num}$ can be defined as a certain percentage of this complete cementation (Supporting Information S3). As an illustration, a simulation case with a cementation of $P_{cem} = 20\%$ is to be interpreted in the sense that, in its initial state, the energy needed to break all its bonds is equal to 20% of that of the fully cemented case.

## 2.3. Granular Fault Gouge Sample

Mineral grains morphologies can be very diverse, but granular gouges generated by comminution are expected to exhibit rather rough and angular shapes (An & Sammis, 1994; Lin, 1999; Olgaard & Brace, 1983). Many studies (Anthony & Marone, 2005; Mair et al., 2002; Nouguier-lehon et al., 2003) have shown that using angular and faceted shapes instead of circular grains led to higher friction coefficients and different shearing behaviours. For 3D laboratory experiments with real grains, steady-state effective friction is usually around 0.6 (Mair et al., 2002), in opposition to spherical particles with effective friction that rarely exceeds 0.45. Because of their invariance by rotation, smooth spherical shapes tend to roll to accommodate deformation of the grain assembly whereas interlocking between angular grains tends to promote dilation. It was also shown recently that mechanical effects of grains surface roughness can only be mimicked by intergranular friction to a certain extent, and that proper modelling of the shear behaviour of granular samples requires shapes with a realistic roughness (Mollon et al., 2020). For these reasons, angular and faceted particles are employed in this study. In order to validate this choice, two samples were implemented: one with angular and faceted particles (used in the remainder of this study) and one with circular grains. This second sample was then discarded, as explained in (Supporting Information S4).

The Matlab package Packing2D is employed to create a granular sample with angular and faceted shapes and the presence of micro-roughness. As a

difference with other studies using multi-sphere approach or clumps of particles (Cho et al., 2007; Potyondy & Cundall, 2004), Packing2D is based on a Fourier–Voronoï method and generates a set of angular and faceted grains with user-defined size distribution and control on key morphological descriptors (such as elongation, circularity, and roundness). This control is performed by choosing a Fourier spectrum that quantifies the frequencies and amplitudes of the grain surface asperities (Mollon & Zhao, 2012). Since the morphological descriptors of the grains of granular gouges may vary significantly between faults, we calibrate our spectrum by visual comparison with published pictures of real gouges (Fig. 3a, b).

The thickness of the shearing zone is ranged between 1 and 5 mm (Rice & Cocco, 2002). We create a 2 mm-thick granular fault gouge (before compaction), resulting in 1.7 mm after compaction, and determine what length is needed to obtain a Representative Surface Element (RSE), (Supporting Information S5). A gouge of $2 \times 20 \text{ mm}^2$ is found to be satisfactory and falls within the same order of magnitude as previous studies (Dorostkar et al., 2017; Ferdowsi, 2014). A fractal size distribution is chosen to comply with the literature on granular gouge composition (Billi, 2005; Billi & Storti, 2004; Blenkinsop, 1991; Muto et al., 2015; Olgaard & Brace, 1983) with a fractal dimension factor D of 2.65. The gouge is composed of 4960 particles with a corresponding equivalent diameter ranging from 28 to 226 µm (average value of 81 µm) and a ($D_{50}$) equal to 70 µm (Fig. 3c).

## 2.4. Numerical Setup for Direct Shear Simulations

Figure 4 presents the DEM model of the simulated granular fault gouge, with rock walls at the top and bottom sides of the granular sample. Contact surfaces of rock walls are sinusoidal to introduce a certain roughness and avoid wall-slip effects since we want to ensure that slip accommodation takes place within the gouge. This is a pure modelling choice, and the rock surfaces should be considered as nominally flat. Inter-particle friction is equal to 1 at the contact interface between walls and particles to make sure that the motion is fully coupled at the wall-
<figure><img src="image_4.png" /><figcaption>Figure 3</figcaption></figure>

a Zoom on a granular sample generated with packing2D, with angular and faceted shapes, and compacted in MELODY2D. b Quartz, photomicrograph (crossed polars) of ATTL fault gouge from (Muto et al., 2015). c Size distribution of the synthetic sample (fractal distribution)

<figure><img src="image_5.png" /><figcaption>Figure 4</figcaption></figure>

DEM model of a granular fault gouge, 4960 angular particles in a 1.7 × 20 mm domain

grains transition. For inter-particle contacts $\mu_{num}$, friction is set to 0.5. This value is often used in DEM gouge experiments and is in the range of inter-particle friction found for two frictional particles in contacts (Kasyap & Senetakis, 2020; Sandeep & Senetakis, 2019). The main parameters used in the numerical model are gathered in Table 1.

The lower wall is fixed, while a normal stress of 40 MPa and a sliding velocity are applied on the upper rock wall. Gravity is ignored in the model, assuming that the fault can be oriented in a wide range of directions, and that gravity forces are negligible compared to those related to normal and deviatoric stresses applied on the gouge. Periodic boundary conditions are present on both right- and left-hand sides of the sample to maintain the continuity of the movement at large slips. The movement

of the upper rock wall in the Y-direction remains free to allow gouge dilatancy.

A dry contact model is considered here, to investigate cementation and porosity influences alone, without fluid. Fault gouge and rock walls are considered rigid with a constant numerical stiffness of 10¹⁵ Pa/m used to limit interpenetration between grains (identical for normal and tangential directions) while mimicking the local deformation of the grains in the contact vicinity (Sects. 2.1 and 2.2). This constant value was chosen to obtain the overall deformability of the sample of the same order of magnitude as the one for bulk granite or shearing modulus (≈ 10–25 MPa depending on initial porosity, see Sect. 3.1). Typical grains interpenetrations under the applied normal stress and with the chosen contact stiffness remain below 1% of the typical grain
<table><caption>Table 1</caption><caption>Numerical parameters and materials properties</caption><thead><tr><th>Property</th><th>Associated variable</th><th>Value</th></tr></thead><tbody><tr><td>Normal stress</td><td>σ<sub>n</sub></td><td>40 MPa</td></tr><tr><td>Shear velocity</td><td>V</td><td>1 m/s</td></tr><tr><td>Rock density</td><td>ρ<sub>r</sub></td><td>2600 grains kg/m<sup>3</sup></td></tr><tr><td>Contact stiffness</td><td>k<sub>n</sub>, k<sub>t</sub></td><td>10<sup>15</sup> Pa/m</td></tr><tr><td>Cohesive bond strength</td><td>C<sub>num</sub></td><td>0 – 2500 MPa</td></tr><tr><td>Percentage of cementation</td><td>P<sub>cem</sub></td><td>0 – 100%</td></tr><tr><td>Percentage of initial porosity</td><td>P<sub>pore</sub></td><td>11% and 16%</td></tr><tr><td>Inter-particular friction</td><td>μ<sub>num</sub></td><td>0.5 (grains – grains) /1 (grains – walls)</td></tr><tr><td>Sample size</td><td>w<sub>gouge</sub> * L<sub>gouge</sub></td><td>1.7 × 20 mm</td></tr><tr><td>Particle equivalent diameter</td><td>φ<sub>min</sub> – φ<sub>max</sub></td><td>28 – 226 μm</td></tr><tr><td>Number of particles</td><td>N</td><td>4960</td></tr><tr><td>DEM time step</td><td>Δt</td><td>10<sup>-9</sup> s</td></tr><tr><td>Proximity updating period</td><td>Δt–contact</td><td>10<sup>-7</sup> s</td></tr></tbody></table>

size. The tangential and normal numerical stiffnesses are equal in this model and grain comminution is disregarded. In this study, we choose to simulate a density of 2600 kg/m³ for particles, leading to an appropriate time step for these simulations of 10⁻⁹ s. An explicit solver is used (Symplectic Euler scheme), to integrate in time the motion of each body. The first simulation stage consists of compacting the gouge to obtain a mechanically stable packing of grains with a controlled granular density. Two different initial packings are considered for this study (i.e. different initial states of porosity $P_{pore}$, (see Supporting Information S6)), mid-dense samples ($P_{pore}$ = 16%) and dense samples ($P_{pore}$ = 11%).

After compaction and stabilization, the contact law between grains is modified by the introduction of cohesive bonds between particles (see the previous subsection) and a 1 m/s slipping velocity is applied on the upper rock wall. This high velocity allows the simulations to run in a reasonable time duration while avoiding disturbing inertial effects since the dimensionless inertial number, in that case, is close to 10⁻³ [i.e. quasi-static dense granular flow according to Da Cruz et al., (2005)].

## 2.5. Numerical Setup for Biaxial Simulations

The strength of granular materials (may they be cemented or not) and rocks is generally considered to follow the Mohr–Coulomb criterion (Handin, 1969; Jaeger, 1971), which states that the maximum tangential stress $τ_{max}$ that the material can withstand, is the sum of a constant term C (generally called “cohesion”) and an additional term tanφ (generally called “internal friction”) multiplied by the normal effective stress:

$$τ_{max} = C + σ_n tanφ \quad (8)$$

Note that the effective stresses are equal to the total stresses as the material is dry. These parameters are not to be confused with the interparticle cohesion ($C_{num}$) and friction ($μ_{num}$) that are introduced in the DEM contact laws, since they represent a collective response of the granular material to shearing.

To characterize our synthetic cemented gouges (in the sens of Mohr–Coulomb), we ran independent simulations of biaxial compression of samples with the same characteristics as those used for our sheared fault gouge (in terms of grain shapes, size distributions, initial solid fraction, interparticle cohesion, and interparticle friction). Rectangular granular samples (4 mm wide and 10 mm high) are placed between four rigid walls (Fig. 5a). The lower wall is fixed in displacement, the upper wall is submitted to a constant downwards velocity $V_y$, and the lateral walls are submitted to a confining pressure $σ_3$, fixed in vertical and rotational motions, and free to move horizontally. For each cementation level and each initial porosity, three tests are performed with confining stresses of 10 MPa, 40 MPa, and 80 MPa. The vertical stress $σ_1$ is monitored during vertical compression (series of pictures with relative damage (Fig. 5a)). Biaxial simulations are only used for the cemented material characterisation Sect. 3.1 in order to obtain the internal friction angle and the cohesion.
<figure><img src="image_6.png" /><figcaption>Figure 5
Illustrative results for a numerical cohesion $C_{num} = 1000\text{MPa}$. a Views of experiment and damage evolution for different percentage of vertical shortening under confining stress of 80 MPa. b Major principal stress as a function of vertical shortening and associated Mohr circles, for three confining stresses (10 MPa, 40 MPa, and 80 MPa) for both dense and mid-dense samples</figcaption></figure>

### 3. Simulation Results

#### 3.1. Cemented Material Characterization Under Biaxial Simulations

Figure 5b provides illustrative results for a numerical cohesion $C_{num} = 1000\text{MPa}$ (corresponding to $P_{cem} = 38\%$ in the dense case and $P_{cem} = 40\%$ in the mid-dense case). Results show that dense samples quickly reach a peak strength before softening towards a plateau value, while mid-dense samples are less stiff and reach the same plateau (approximately at the same level of vertical strain) without passing by a well-defined peak. When plotting the Mohr circles corresponding to the maximum values

of $\sigma_1$ (Fig. 5b), we confirm that all samples follow the Mohr–Coulomb failure criterion, which makes it possible to characterize their cohesion $C$ and internal friction angle $\phi$ (in the sense of Mohr–Coulomb). For each density, elasticity moduli $E$ (not strictly equal to Young’s moduli because of biaxial conditions) are also extracted based on the vertical stiffness of the samples at low ($< 1\%$) vertical strains.

Similar Mohr–Coulomb graphs were plotted for the other $C_{num}$ values and the same conclusion is drawn concerning the Mohr–Coulomb criterion. Results are gathered in Fig. 6 and show that internal friction $\phi$ is only moderately influenced by the cementation level in the mid-dense samples, and is a
<figure><img src="image_7.png" /><figcaption>Figure 6
a Internal friction angle (°), b cohesion (MPa) and c elasticity moduli (GPa) as a function of the cementation (%) for dense and mid dense samples</figcaption></figure>

non-monotonous function of the cementation in the
case of dense samples. It is also strongly linked to
initial density: mid-dense samples have friction angle
values ranging from 32° to 35°, while dense samples
have friction angle values ranging from 36° to 44°.
Cohesion increases monotonously with cementation
[as found by Wissa (1965) and (Lade and Overton
(1989)] and reaches much higher values for dense
samples (up to 70 MPa for Pcem ≈ 100%) than for
mid-dense samples (up to 20 MPa). It can also be
noted that a certain percentage of cementation is
needed to have a measurable cohesion within the
gouge (Pcem > 20% for mid-dense samples and
Pcem > 10% for dense samples). However, Elasticity
moduli are rather unaffected by cementation and is
close to 16 GPa for mid-dense samples and to 25 GPa
for dense samples.

## 3.2. Influence of the Cementation on Gouge Kinematics

### 3.2.1 Effective Friction and Dilation

For dense samples, several simulations with various
cementation levels are performed with the numerical
setup for a direct-shear experiment (Sect. 2.4). In the
following sections, we will use the effective friction
μ* = τ/σn (i.e. total tangential force τ resisting to

fault sliding, divided by the applied normal force σ_n)
and the dilation ε_y = dy/w_gouge (i.e. volume variation
dy in the direction of the gouge thickness w_gouge) to
analyze the mechanical behaviour of gouges, as done
in previous studies (GRD Azéma & Radjaï, 2014;
Berger et al., 2015; Midi, 2004). It is important to
keep in mind that cementation makes μ* dependant
on the normal load and thus different from a
Coulomb-like friction coefficient, hence the term
“effective friction”.

Figure 7a, b provides typical curves of the
measured effective friction of the fault as a function
of the horizontal displacement imposed to the upper
wall. In all cases, the tangential force increases
linearly until a maximum effective friction μ*,
demonstrating the maximal effort that the loading
system must provide to overcome interlocking in the
gouge and accommodate imposed shearing. The
linear elastic part of the curves represents the
stiffness of the granular material. It should be noted
that this stiffness is only related to the gouge layer
itself since no other compliance (related to the
surrounding medium, for example) is considered in
the simulations. All simulations with initial dense
samples follow the same elastic slope, and adding
more cementation extends the elastic part before the
peak (Fig. 7b).
<figure><img src="image_8.png" /><figcaption>Figure 7</figcaption></figure>

a Effective friction curve for different percentages of cementation for dense samples as a function of the slip distance (µm). b Zoom in on the peak. c Brittleness as a function of cementation, the brittleness is expressed as the ratio of $(\mu_p^* - \mu_{ss}^*)$ divided by $\mu_p^*$. d Dilation variation (µm) as a function of the slip distance (µm) for dense samples. e Dilation variation (µm) as a function of the slip distance (µm) for dense samples, longer slip distance. Total simulation from "A" to "H". Letters correspond to different steps in curves presented: [A] is the initial state before shearing and is identical in all cases. [B] is the peak location, [C] is the end of the first peak (only appears for dense samples with a lot of cementation), [D] is half-peak, [E] is the end of the major dilation phase, [F] is the observed end of the effective friction peak, [G] is a common state for beginning of steady-state and [H] is the end of simulation with s = 2500 µm (results table in Supporting Information S7). The case "0 cementation" is not represented here as it is similar to the case with 4% cementation in terms of friction and dilation variation. It should be noted that these results are not in the range of small deformations, as the total slip displacement (H = 2.5 mm) is higher than the gouge thickness (≈ 1.75 mm)

According to the literature, shear strength is enhanced with cohesion, and thus with the increase of cementation (Wissa, 1965), (cf. Fig. 6b).

Depending on the cementation percentage of each simulation, the strength peak evolves from a smooth, delayed, and of moderate amplitude (poorly
cemented cases) to a sharp, short, and higher
amplitude (highly-cemented cases), (Fig. 7a, b).
Cohesive strength in granular materials is known to
correlate with an increase in the brittleness of the
material (Das et al., 2014), and this is confirmed by
(Fig. 7c) which presents the evolution of brittleness
within the sample (Meng et al., 2021). A qualitative
change in behaviour is observed close to
$P_{cem} = 10\%$. This transition corresponds to the value
where cementation starts to induce cohesion within
the gouge (red line). From this percentage and above,
a marked evolution in terms of peak strength is
observed, with the emergence of a sharp peak
strength increasing with cementation. After the peak,
the effective friction decreases in all cases towards a
plateau, and does not evolve significantly until the
end of the simulation. All steady-state values oscillate
around the same effective friction $\mu_{SS}^*$, (averaged
from the beginning of the plateau until the end of the
simulation), close to 0.5 (ranging from 0.45 to 0.51).
It is interesting to notice that this value is in
agreement with other numerical studies (Rathbun
et al., 2013), but lower than typical 3D experimental
values (which are usually above 0.6). This discrep-
ancy is related to the 2D character of the simulations
(Frye & Marone, 2003).

For dilation, different behaviours are also
observed according to the percentage of cementation
within the gouge (Fig. 7d). For $P_{cem} \le 10\%$, the
behaviour is similar to the case without cementation
as the cohesion in the sense of Mohr-Coulomb is
smaller than 1 MPa (which is negligible compared to
the normal stress). A progressive dilation is observed
until the steady-state “G”. The increase of cementa-
tion (for $P_{cem} \ge 10\%$) seems to accelerate the dilation
of the gouge (Fig. 7d), which appears earlier with
regard to the slip distance. This dilation enhancement
is possible because some frictional contacts emerge
with the breakage of cemented bonds. Observing the
second part of simulations from $s = 1$ mm to
$s = 2.5$ mm (Fig. 7e), three main tendencies are
highlighted: materials with medium cementation
($\sim 10\% \le P_{cem} \le 75\%$) reach a stabilized density
whereas poorly cemented materials continue to dilate
and highly cemented materials present a contractive
behaviour. Poorly cemented materials ($P_{cem} \le 10\%$)
present more frictional contacts than cohesive

contacts and frictional particles are known to enhance
dilation (Roy & Luding, 2017). For high cementation
($P_{cem} \ge 75\%$), although the steady-state friction value
announces a friction stabilization “G”, a contracting
behavior is observed later (from “G” to “H”). This
contracting is attributed to the rupture of the asper-
ities formed by cohesive links (i.e. agglomerates) at
the beginning of the shear localization “F”, (details
are given in Sect. 3.4).

3.2.2 *Interface Failure Modes*

Cohesive strength (correlated with brittleness
increase) has a major influence on strain localization
(Maurer, 1965). From observations on friction and
dilation (for dense samples cases), our results can be
gathered in three different cemented materials: poorly
cemented ($P_{cem} \le 10\%$), cemented
($10\% \le P_{cem} \le 75\%$), and highly-cemented
($P_{cem} \ge 75\%$), corresponding to three modes of
deformation zone. Figures 8, 9, and 10 present the
evolution of damage between grains during sliding
for each of these three states of cementation. The
rupture of a cohesive bond is represented by the
increase of damage to the concerned grains. This
damage is set to 0 for each grain when cohesive
bonds are first established (all the bonds are intact)
and may evolve until 1 if all these bonds reach the
“broken” status (cf. Sect. 2.2). It is thus a relative
damage with respect to an initial state. The repre-
sentation of the relative damage gives a picture of the
state of cementation between grains and their location
within the gouge. As shear band formation and
evolution varies with the level of bond strength (Jiang
et al., 2013), the damage will allow following the
formation of failure patterns (Riedel cracks, shear
bands, etc.) and their orientation.

(i) The poorly-cemented material is a medium where inter-particular cohesion introduced between grains is not sufficient to maintain cohesive bonds during shearing. Figure 8 presents the simulation for $P_{cem} = 4\%$, revealing very few tensile forces in accordance with the low percentage of cementation in the sample. As soon as the upper rock wall is set into motion, almost all cohesive bonds break and only a few of them resist until friction
<figure><img src="image_9.png" /><figcaption>Figure 8</figcaption></figure>

Relative damage snapshot for $P_{cem} = 4\%$ in the dense sample (entire granular gouge). Letters correspond to different steps in the curves presented in Fig. 7 (A4: no damage and no-slip distance, G4: maximum damage and slip distance of 1 mm). (i) Tensile contact force orientations allow getting a grasp on the way load is transferred through a granular sample. This information is coded in a polar diagram which provides the distribution of the orientation of tensile contacts

peak “B₄”. Dilation is found to influence shear bands and strain localization (Mead, 1925), and the limited dilation in the first stages of shearing is consistent with the absence of shear band formation.

(ii) The second type of behaviour in Fig. 9 corresponds to a cemented granular material ($P_{cem} = 38\%$). In contrast with the first case, this material presents clear augmentation of tensile forces for both initial states “A₃₈” in dark blue and final states in red in Fig. 9a, b. Tensile contacts then reduce from effective friction peak “B₃₈” to the end of the peak “F₃₈” with the breakage of cohesive bonds, but many of them remain intact, and influence gouge behaviour. Friction peak “B₃₈” is the starting point of a movement within the gouge with the highlighting of a preferential localization of cohesive bonds rupture [white arrow in Fig. 9 is considered as the first Riedel deformation (Tchalenko, 1970)]. This rupture develops in the next stages of the

simulation “B₃₈ to H₃₈” with a pattern similar to a Riedel crack R (oriented in the sheared direction, ~ 12° from the upper wall). The progression of the Riedel cracks towards a shear band increases until the end of the effective friction peak, where it is no longer detectable among the damage zone. The different Riedel geometries are associated with different shear deformation degrees. R' shear bands are not visible in the numerical results, but we observe tensile fractures T. This is in line with findings reported in Cho et al., (2008), who highlighted that tensile fracturing is very important in the development of a shear zone, and particularly at low normal stress. From friction peak “B₃₈”, damage evolution also highlights the presence of cohesive agglomerates formed by intact bonds within the gouge (ex: Ag₁), with size decreasing with time. In contrast to poorly cemented material where contacts take place between two particles, contacts occur here between big
<figure><img src="image_10.png" /><figcaption>Figure 9</figcaption></figure>

Relative damage snapshot for $P_{cem} = 38\%$ in the dense sample (entire granular gouge). Letters correspond to different steps in the curves presented in Fig. 7. (A38: no damage and no slip distance, H38: slip distance of 2.5 mm). “Ag1” denotes a large agglomerate followed during its motion. The white arrow follows the low angle Riedel shear “R” inside the gouge and “T” is an example of extension (or tensile) fracture.

(i) Tensile contact forces orientation allows getting a grasp on the way load is transferred through a granular sample. This information is coded in a polar diagram which provides the distribution of the orientation of tensile contacts

agglomerates of cohesive grains. These agglomerates can lead to rather inhomogeneous behaviours inside the granular gouge, changing the whole geometry and particle size distribution. Similar observations were made by (Cho et al., 2008) and Rognon et al., (2008). Cohesive agglomerates also participate in increasing the overall dilation before steady-state (Fig. 7) as found by Lade and Overton (1989) for large cemented particles at low confining pressures. Even though dilation is needed to observe shear bands, this cemented material (presenting the highest dilation after peak strength) inhibits the

persistence of shear bands until the end of the shearing. This may be due to the important breakage of cohesive links into frictional contacts after the end of effective friction peak.

(iii) Increasing again $P_{cem}$ leads to a highly-cemented zone ($P_{cem} = 95\%$) where most cohesive bonds stay intact during the entire simulation (Fig. 10). The numbers of tensile forces are obviously higher than in previous materials (Fig. 10i) and correspond to Lade and Overton (1989) results. This highly cemented rock enhances the formation of two Riedel cracks R at the second friction peak in “$C_{95}$”.
<figure><img src="image_11.png" /><figcaption>Figure 10</figcaption></figure>

Relative damage snapshot for $P_{cem} = 95\%$ in the dense sample (entire granular gouge). Letters correspond to different steps in the curves presented in (Fig. 7). (A95: no damage and no slip distance, H95: slip distance of 2.5 mm). R represents the Riedel shear bands and Y the horizontal shear localization. (i) Tensile contact force orientations allow getting a grasp on the way load is transferred through a granular sample. This information is coded in a polar diagram which provides the distribution of the orientation of tensile contacts

Then, “D95” shows the largest Riedel band-
width, as this step also corresponds to the
maximum dilation peak in the simulation. In
the next steps, from “E95” to “H95”, Riedel
crack thicknesses progressively reduce in favour
of a horizontal shear localization S at the top of
the granular gouge in “H95”. This sheared
thickness reduction was also observed in other
numerical studies with the increase of strain
inside the model (Cundall, 1989). As previously

supposed (Sect. 3.2.1), the progressive breakage
of cohesive agglomerates (i.e. asperities) at the
interface between gouge and rock boundary
seems to explain the contractive response
observed in (Fig. 7e) from “G95” to “H95”.
The shearing localizes on the bottom or top part
of the gouge (depending on the simulation) and
lets the majority of particles (still cemented)
behave as a single solid-like body. This shearing
or “active” zone increases with slip distance.
### 3.3. Influence of Initial Porosity with Cementation

The same percentages of cementation have been tested on a second set of samples with a larger initial porosity $P_{pore} = 16\%$ (termed as mid-dense in what follows, as opposed to the dense samples with $P_{pore} = 11\%$). The three cementation levels previously highlighted have been kept to compare the results between dense and mid-dense samples. The trend observed for mid-dense samples is similar to dense samples, but an influence of both initial porosity and cementation on the strength of cemented materials is noted.

The increase of initial porosity is supposed to reduce the fracture strength (Taylor, 1948). Our results present the same conclusion, but different shapes of effective friction peak are observed: the peak strength may be sharp, short, and intense (dense and highly cemented cases) or smooth, delayed and of moderate amplitude (mid-dense and moderately cemented cases), (Fig. 11a, b). The denser sample shows higher dilation rates than the mid-dense sample (Fig. 11d), as the initial gouge is initially more compacted (Wood, 1990). When the steady-state is reached, maximum dilatancy is also observed for denser initial samples with twice the deformation obtained with mid-dense samples, for all initial cementations ($\epsilon_{y-SS} = 8\%$ for dense samples vs $\epsilon_{y-SS} = 4\%$ for mid-dense).

Figure 11c gathers the brittleness based on strength ratio (Bishop, 1971) as a function of cementation for the whole simulation campaign and gives more precise information on the transition zone. For samples where a measurable Coulomb cohesion exists, a higher percentage of cementation increases the brittleness of the sample (Fig. 11c) for both dense and mid-dense samples, and thus leads to enhance the strength of the gouge (Fig. 11a). However, the increase of initial porosity leads to a weaker material in terms of cohesive strength. Without the effect of cementation, the brittleness is almost four times higher for dense samples (i.e. ~ 0.4) than for mid-dense samples (i.e. ~ 0.1). The mid-dense case with the lower cementation is, here, the weaker material tested. When the cementation level reaches values close to its maximum ($P_{cem} \approx 100\%$), then the brittleness can reach almost 0.9 in the dense material

and almost 0.66 in the mid-dense material. It is very likely that such high values are related to the fact that the cementation of the sample is very high while the confining stress (40 MPa) is rather low (Lade & Overton, 1989).

Mid-dense samples also present higher relative damage with easier grains reorganizations (Fig. 12), which agrees with the relatively low Coulomb cohesion involved. One high angle Riedel shear band $R'$ in the 40% case, and two high angle Riedels $R'$ in the 95%, are followed by a horizontal shear localization $Y$ at the bottom of the granular gouge Riedel bands observed do not have the same orientation as those observed with dense materials and are less persistent during shearing. The more ductile character obtained for mid-dense samples could be a reason for the different Riedel angles observed (Misra et al., 2009). The non-persistence of shear bands for mid-dense samples is in adequation with previous studies on sandstones (Antonellini & Pollard, 1995; Dunn et al., 1973). They have shown that gouges with high porosities enhance a distributed deformation that takes the form of cataclastic flow, in contrast with low porosity sandstones which better fail by localization and strain softening (Hirth & Tullis, 1989).

### 3.4. Influence of Ductility with Cementation and Porosity

Force chains are key elements to understand the kinematics of the gouge layer and illustrate changes in ductility behaviour inside the granular gouge as a function of cementation or porosity. It should be noted that, in the present context, “ductility” is to be understood in the sense of granular mechanics: it is related to the suddenness or slowness of the post-peak frictional weakening but does not involve any visco-plastic phenomenon.

#### 3.4.1 Evolution of Ductility with Cementation

Commonly used in granular physics (Majmudar & Behringer, 2005; Zhang et al., 2017) for load transfer observation, force chains are plotted between grains centres (for each contact). The norm of the contact force is coded by the thickness and colour of each chain (in case of edge-edge contacts, handled through
<figure><img src="image_12.png" /><figcaption>Figure 11</figcaption></figure>

a Effective friction curve, comparison between dense and mid-dense samples as a function of the slip distance (µm) for the three different cemented material highlighted. b Zoom in on the peak for mid-dense samples. c Brittleness as a function of cementation for both porosity states, the brittleness is expressed as the ratio of $(\mu_p^* - \mu_{ss}^*)$ divided by $\mu_p^*$. d Dilation variation (µm) as a function of the slip distance (µm) for dense and mid-dense samples. Letters correspond to different steps in curves presented: [A] is the initial state before shearing and is identical in all cases. [B] is the peak location, [C] is the end of the first peak (only appears for dense samples with a lot of cementation), [D] is half-peak, [E] is the end of the major dilation phase, [F] is the observed end of the effective friction peak, [G] is a common state for beginning of steady-state and [H] is the end of the simulation (results table in Supporting Information S7)

the nodal contacts at the extremities of the contact segment, the total resulting force is considered and coded in a unique force chain). It allows to get a grasp on the way load is transferred through a granular sample. Similar information can be coded in polar diagrams which provide the distribution of the orientation of contact normal vectors. The poorly cemented materials ($P_{cem} < 10 - 25\%$) show a ductile behavior, with a progressive particle reorganization, and no localized shear is observed in Fig. 8. It

displays a preferential orientation for force chains inclined at 45° from the upper rock wall (Fig. 13a). Similar orientations of force chains have been observed by Morgan and Boettcher (1999) for circular particles. The change of orientation of normal forces at 45° appears before the effective friction peak, showing that the gouge started to dilate before reaching the peak. The evolution of the granular flow gives way to simple contact law with inter-particle friction only. Once the gouge has
<figure><img src="image_13.png" /><figcaption>Figure 12</figcaption></figure>

a Relative damage snapshots for mid-dense sample (entire granular gouge). For both peak [B] and end [F] of effective friction peak. Letters correspond to different steps in the curves presented in Fig. 11. R' represents high angle Riedel shear bands and Y the horizontal shear localization. Evolution of ductility with cementation and porosity. b Tensile contact forces orientation allows to get a grasp on the way load is transferred through a granular sample. This information is coded in a polar diagram which provides the distribution of the orientation of tensile contacts

dilated, grains can reorganize to allow shearing, and the gouge tends towards a steady state of sliding from the end of the effective friction peak to the end of the simulation. The combination of limited cementation (mostly frictional contact, no cohesion in the sense of Mohr–Coulomb), gouge dilation (Fig. 7), and preferred orientation of force chains confirm a typical granular Couette flow for the poorly cemented zone (GRD Da Cruz et al., 2005; Midi, 2004).

With the increase of cementation ($P_{cem} = 38\%$), the orientation at 45° is not yet present at friction peak B, because no significant motion has occurred yet and grains are still in their initial configuration. Instead, there is a large number of smaller ramified force chains distributed in a homogenous distribution of force networks. The strength of the material is linked to the development of force chains (Zhang et al., 2017), and our results confirm that the intensity of force chains is enhanced with cementation, as
<figure><img src="image_14.png" /><figcaption>Figure 13</figcaption></figure>

a Snapshot of the force-chains network for dense samples at friction peak [B] and the end of simulation [H] (from 0 to $10^4$ N). Force chains are defined according to the typical level of stress present in the system (here the applied normal stress) and thickness evaluates between 0 and 1 between the minimum and maximum strength (0–1e4 N). A cut-off is made at 1e4 N to have the same visualization for all cases (higher force chains are observed for a high percentage of cementation). b Zoom in on the material with $P_{cem} = 95\%$, snapshot of the force-chains network at the end of friction peak [F] and the end of the simulation [H] (from 0 to $10^4$ N). (i) Tensile contact forces orientation allows getting a grasp on the way load is transferred through a granular sample. This information is coded in a polar diagram which provides the distribution of the orientation of tensile contacts

strength is also enhanced (Fig. 13a). Combining force chains and relative damage makes it possible to identify from the force network the cemented agglomerates highlighted in Sect. 3.2.2 (red arrows). Each cemented particle is crossed by a force chain, while isolated frictional particles are not (Fig. 13a).

These agglomerates also modify the particle size distribution within the gouge, which is well-known to act on shear bands formation in addition to the initial density of the sample (Marone & Scholz, 1989).

Figure 13a, b displays a dense sample with 95% cementation at the end of friction peak. Some
particles are detached from the cohesive band and operate alone (within the active or shearing zone), defining clear geometrical asperities. These asperities create contacts between the cemented layer and the rock wall, which concentrate the normal load. A steady-state is reached when enough of this tribological third body has been released to avoid any asperity contact and to produce a three-body sliding. The third body considered here is the gouge material within the shearing contact zone. More information on the third body concept is available in Iordanoff et al., (2005) and Fillot et al., (2007). Gouge first dilates thanks to the emergence of asperities (i.e. grains emerging from cohesive bond breakage) and agglomerates of grains, forming inter-particle bridges, which releases gouge in the interface as the sliding progresses. Normal forces orientation follows the same trend as the dilation: a first increase at the first step (due to inter-particle bridge formation) and a decrease back to the initial state, as the majority of cohesive bonds are still active (Fig. 13i). Cementation here also enhances local dilation during the strength weakening, increasing local porosity inside the gouge as found by Faqih et al., (2006).

### 3.4.2 Evolution of Force Chains with Initial Porosity

Figure 14 compares force chains and normal contact orientation for dense and mid-dense samples. Whatever the cementation level, the initial preferred orientation of normal forces at zero displacement (step A) changes between dense and mid-dense samples (Fig. 14b). Dense samples lead to a homogeneous repartition of contact normal orientations and mid-dense samples present normal forces mostly oriented perpendicularly to rock walls. The evolution between the three cemented materials at effective friction peak is similar for both values of porosity, but samples initially denser favour a more homogenous distribution of the force network [i.e. force chains (Fig. 14a)] passing through almost all particles (consistent with previous section showing that the strength increased with a highly cemented and dense sample). We can assume that this difference is related to the more brittle behaviour of the dense sample, for which no reorganization and dilation of grains occurs before the peak. It preserves the initial homogeneous

contact network, but makes the failure more sudden, in contrast with the mid-dense case where irreversible grains motions occur before the peak.

## 4. Discussion: A Link Between Gouge Rheology and Fault Zone Structure

### 4.1. Shear Bands and Weakening Mechanisms

#### 4.1.1 Shear Bands and Critical Dilation

Critical dilation can be defined as the necessary dilation for macroscopic shear failure formation (Kranz & Scholz, 1977). In our simulations, it corresponds to the dilation $ε_{yp}$ obtained at the effective friction peak (Figs. 7, 16) which is different from the final dilation $ε_{y-SS}$ previously introduced.

* For cases with limited cementation ($P_{cem} < 10\%_{dense} - 25\%_{middense}$) — where 100% frictional contact is quickly reached and no localized shear is observed — the critical dilation is similar for both dense and mid-dense samples and presents the highest values (~ 2% in Fig. 15). A high dilation rate is present in the pre-peak phase and explains the important dilation observed at the friction peak. Porosity alone (without cementation effect) does not have a major influence on critical dilation.

* For cases with moderate to high cementation ($P_{cem} > 10\%_{dense} - 25\%_{middense}$), two main trends are observed with the increase of cementation: (1) a decrease of the critical dilation for mid-dense samples, (2) an enhancement of the critical dilation for dense samples. It can be noted that all critical dilations observed in this zone are smaller than in the rectangle “i”. In cohesionless materials, dilation is supposed to be larger for a denser granular gouge, but the introduction of cementation creates the opposite effect. Similar results have been found for other dense cohesive materials through other types of experiments (Faqih et al., 2006). The critical dilation is here smaller for dense samples than for mid-dense samples, a result that may be due to the real cohesion values involved within the sample (i.e. higher for dense sample).
<figure><img src="image_15.png" /><figcaption>Figure 14</figcaption></figure>

a Snapshots of the force-chains networks at the effective friction peaks [B] for six chosen cases (0–1e4 N). For 3–4, 40–38, and 95% of cementation. In the pictures, only a quarter of the total gouge is displayed, and even if the global behaviour is similar to one quarter, force chains are not homogeneously distributed within the gouge after friction peak. b The graph provides the numbers of contacts with a normal vector oriented in a given direction (using polar diagrams where the angle-axis is the orientation and the radius-axis is the number of contacts), and with a colour evolving during the different stages of each simulation

* Comparing shear band formation and critical dilation, we can deduce that with more than 1% of critical dilation (for both dense and mid-dense samples), shear band formation is inhibited. Then, shear bands persistence until steady-state depends mostly on cementation, and occurs for highly-cemented materials, corresponding to a critical dilation between 0.08% and 1% (rectangle “iii”). (Noel et al., 2021) found that for Tavel limestone (11% initial porosity corresponding to our dense sample), a critical dilation of about 1% has to be reached to obtain shear failure of an initially ductile rock subjected to a pore fluid pressure

increase. Their higher critical dilation could be explained by the fact that they study intact rocks, but also by the higher confining pressure used and the presence of fluid.

### 4.1.2 Evolution of Weakening Mechanisms

In Fig. 16, peak length $D_{pp}$ is the distance between the peak of effective friction and the beginning of the steady-state. $D_{pp}$ value is relevant for rupture energy calculation and can be related to the energy needed to weaken the fault until the steady-state of sliding:
<figure><img src="image_16.png" /><figcaption>(b) Critical dilation</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15</figcaption></figure>

a Total dilation as a function of slip distance, the marked point locates the critical dilation for the 3 cemented materials. b Critical dilation for macroscopic failure for both dense and mid-dense samples as a function of the percentage of cementation

* (Fig. 16c) highlights two main trends for $D_{pp}$: (i) with $P_{cem} < 10\%_{dense} - 25\%_{middense}$, mid-dense samples need a higher displacement than dense samples to reach their steady-state regime (810 vs 300 µm), (ii) with $P_{cem} > 10\%_{dense} - 25\%_{middense}$, $D_{pp}$ decreases with $P_{cem}$, similarly for both dense and mid-dense samples. Initial porosity seems to affect the peak duration, but only for cases with limited cementation. When cementation is large enough, the initial state of porosity doesn't have an additional influence on $D_{pp}$. The transition between the two zones corresponds to the mechanical limit between the cohesionless and cohesive model (i.e. formation of cohesion in the sense of Mohr-Coulomb theory) represented by the red line on graphs. Adding more cementation within the gouge seems to diminish the porosity effect regarding the peak length. The evolution of the peak distance $D_p$ (i.e. distance reached at effective friction peak), is similar to the evolution of critical

dilation (Fig. 16b), and increasing the cementation also tends to balance the change in porosity, leading to a similar peak location.

* Double weakening shapes are observed in some effective friction curves in Fig. 16. This pattern seems to be a characteristic of dense cohesive material as our mid-dense case only has one weakening phase. When it exists, the first weakening period seems to follow a decreasing affine law and appears for dense samples with $P_{cem} \ge 10\%$. Isolating contributions from each mechanism (bond breakage, dilation, friction), it appears that this first weakening period is mainly due to interparticle friction and cement breakage, generating high frictional contacts between the few grains which are not cohesive anymore. The second weakening period obtained for dense cemented material (and corresponding to a simple weakening for mid-dense material), is mainly due to dilation mechanisms. This energetic separation in two phases adds a new key for understanding weakening mechanisms. Similar double weakening has been observed experimentally as reported in the paper from (Paglialunga et al., 2021).

The weakening zone usually appears with the formation of B and R shears bands (Mair & Marone, 1999). However, with the separation into two weakening phases, it appears that Riedel bands only form at the beginning of the 2nd weakening phase (letter C in previous figures). In fact, the first weakening phase just allows movement within the gouge by breaking some isolated bonds, but actual sliding needs dilation to operate and thus is only observed within the 2nd weakening.

## 4.2. Cemented Materials and Mohr-Coulomb Theory

### 4.2.1 Comparison to Real Cemented Materials and Rocks

The cementation considered in this paper has a clear definition from a numerical point of view but is difficult to characterize on real samples. For this reason, it is instructive to consider instead the Mohr-Coulomb properties of the cemented gouges, for which characterization techniques are available
<figure><img src="image_18.png" /><figcaption>Figure 16</figcaption></figure>

a Effective friction curves, b peak location $D_p$ (μm), c peak length $D_{pp}$ (μm), for both dense and mid-dense samples, as a function of the percentage of cementation within the gouge

(Sect. 2.5). It can be noted that the samples used are large enough to avoid any size effect or any thickness dependence. Figure 17 compares the different

synthetic materials tested in this paper (with different percentages of cementation and porosities) with typical rock properties found in the literature

<figure><img src="image_19.png" /><figcaption>Figure 17</figcaption></figure>

Internal friction angle as a function of Mohr–Coulomb cohesion for the simulated materials, as characterized by numerical biaxial tests (Fig. 5). Dense samples are marked by squares and mid-dense samples by diamonds. The discs correspond to experimental characterizations of different rocks referenced in (Schellart, 2000) listed in the table next to the graph
(Schellart, 2000), in terms of internal friction angle and cohesion. Our samples provide adequate values of internal friction angle (especially in the mid-dense case) and are in the lower range of cohesions. The low cohesion is consistent with the fact that we are simulating cemented gouges instead of intact rock.

Figure 17 confirms that our samples have similar Mohr-Coulomb properties to some real rocks: (a) quartz sand with poorly cemented dense samples, (b) rock “8” (Tennessee sandstone) with 75% cemented gouge (dense), (c) rock “1” (granite) with highly cemented gouge (mid-dense). Even though some friction angle and cohesion are close to some Mohr-Coulomb values, a fault gouge develops a complex Riedel band structure reflecting a heterogeneous stress field different from the one observed for intact rocks (Lockner & Beeler, 2002).

Future work might be undertaken to target the Mohr-Coulomb properties of a cohesive rock for which direct shear experimental results are available, in order to propose a direct comparison with our simulations.

### 4.2.2 Discussion About Mohr-Coulomb Theory

Cemented gouges have interesting behaviour because they are not intact rock and neither are simple frictional gouges. In this section, we compare failure criteria applied to intact rock or fault gouge layer to the cemented gouges presented previously in order to evaluate which criterion can best fit the failure of cemented gouges.

Mohr-Coulomb theory is often used to predict failure of intact rocks (Handin, 1969; Jaeger, 1971). Byerlee and Savage (1992) also analysed a fault gouge layer under shearing and observed that it behaves as a Coulomb material, with an elastic deformation until failure. However, Marone et al., (1992) proposed another adapted criterion that is more suitable to fault gouge layer. The two failure criteria considered to support these assumptions are:

$$ \tau/\sigma_N = C/\sigma_N + \tan\varphi \quad (9) $$

$$ \tau/\sigma_N = C\cos\Phi/\sigma_N + \sin\Phi \quad (10) $$

The first criterion applies to intact rock and the normal and tangential stresses $\sigma_N$ and $\tau$ apply on a

shear crack with a clear orientation. It is called “Coulomb failure” and provides the “true” Mohr-Coulomb properties $C$ and $\varphi$ of the intact rock. The second criterion applies to a domain where shearing occurs in a distributed and continuous way and the stresses $\sigma_N$ and $\tau$ are not requested to apply on a specific crack. Rather, their orientation is related to that of the principal stresses that develop in the sheared domain. This is known as “Coulomb plasticity”. As explained in (Marone, 1995), the first criterion could be applied to faults with very thin layers of gouge, while the second one might be relevant to the case of thick layers of granular gouge where such rotation of the principal stresses could occur. In that case, Coulomb plasticity could lead to a reduction of the effective friction of the fault. This assumption was validated experimentally in a cohesionless case (Marone, 1995).

It is instructive to determine to which extent these two expressions can predict the fault strength, solely from the “true” Mohr-Coulomb properties of the rock. This comparison can be done using our numerical results on fault direct shear tests and biaxial tests. Figure 18 thus presents the effective friction in the fault as obtained numerically with direct shear simulations in Sect. 3 of the present work, as a function of the effective friction predicted by both failure criteria 9 and 10 based on our independent determination of $C$ and $\varphi$. With this graph, we can directly compare the peak of strength ratio obtained with both models (direct shear and biaxial test) realized with the same conditions and characteristics (cf. Sect. 3.1). Indeed, gouge behaving as granular material (i.e. with no cohesion) present lower friction than for intact rock as predicted by literature (points under the dotted line). However, cemented gouges present the opposite behaviour once the percentage of cementation account for a small cohesion (i.e. for more than $10\% < P_{cem}$ for dense samples and $25\% < P_{cem}$ for mid-dense samples). The increase of cementation increases the effective friction of the cemented gouge compared to intact rock, and the reduction of porosity overstate the difference. It appears that the proposed adapted law 10 for fault gouges (Marone et al., 1992), does not fit very well with cemented gouges either. The initial stress ratio from the Coulomb criterion (Eq. 9) is closer to values
<figure><img src="image_20.png" /><figcaption>Figure 18</figcaption></figure>

a Effective friction as measured with direct shear experiments (simulation campaign) as a function of effective friction as predicted by failure criterion (calculated with internal friction angle and cohesion found with biaxial tests), for dense and mid-dense samples. The blue line (empty marker points) represents the effective friction given for sheared gouge with Eq. 10 (Marone et al., 1992) compared to the grey line (filled marker points), effective friction for intact rock with Eq. 9 (Handin, 1969). b Zoom in on the less cemented materials

found for our cemented gouges but is not suitable as cementation increases.

We hypothesize that one of the reasons that could explain these differences could be the different boundary conditions used for the two kinds of tests (direct shear/biaxial experiment). In particular, the presence of the two walls of the fault might act as barriers to the natural development of Riedel shear bands and increase the fault strength. Another key of understanding could be the heterogeneity of the stress field in the fault (in relation to the kinematic constraints mentioned above), which would prevent a simplistic use of the criterions 9 and 10. As we need further analyses and simulations to confirm these hypotheses, we shall discuss these propositions in future work comparing different gouges with cement or matrix.

## 5. Conclusion

A 2D granular fault gouge model has been implemented in the framework of DEM in order to establish a link between gouge properties (cement and initial porosity) and rheological behaviour. Our results confirm that both initial porosity and cementation influence gouge behaviour. They play a role in the internal structure and geometry of the gouge and thus modify the rheological behaviour of the fault. In the range of our numerical experiments (normal stress 40 MPa and slip velocity 1 m/s), the increase of cementation within the gouge leads to an increase of effective strength whereas an increase of initial porosity tends to a reduced strength.

(i) The strength peak evolves from a smooth, delayed, and of moderate amplitude (mid-dense and poorly cemented cases) to a sharp, short, and higher amplitude (dense and highly-cemented cases). Brittleness is enhanced with cementation, especially in the case of dense materials. It
evolves similarly with the internal cohesion within the gouge and leads to different failure patterns. The highly-cemented material with low initial porosity presents clear Riedel band formation evolving as in the theory of Tchalenko (1970). The high cohesion and internal friction angle values make this material very close to intact rock properties. In contrast, the weakest material is obtained with the highest porosity state and for no cementation: no failure pattern is observed and a typical granular Couette flow is highlighted.

(ii) Shear failures need a critical dilation to form. For poorly-cemented material, not presenting enough cohesion (i.e. $C < 1$ MPa for poorly cemented material), the same critical dilation is observed for both dense and mid-dense samples. The initial porosity does not have a major influence on critical dilation for cohesionless materials. However, when cohesion is present within the particles, the critical dilation appears to be smaller for denser samples and increases with cementation, whereas mid-dense samples present the opposite behaviour.

(iii) Effective friction curves present double weakening shapes for dense samples with enough cementation. The first weakening phase is triggered by the rupture of cohesive bonds and friction, whereas the second weakening comes with the dilation peak. The peak length (analogous to the $D_c$ from rate and state laws) mainly decreases with the percentage of cementation.

(iv) The increase of cementation increases the effective friction of the cemented gouge compared to intact rock, especially with dense cemented gouges. Effective friction, from the Coulomb failure criterion or the Coulomb plasticity criterion (Marone et al., 1992), does not predict very well our values for cemented fault gouges. The different boundary conditions and stress fields obtained for the direct shear experiment and the biaxial test could explain the difference, but further analyses are needed to confirm these hypotheses.

In this paper, we simulated a cemented granular fault zone with angular and faceted grains linked by

cohesive bonds to model the cementation. Although we can easily measure the percentage of cementation within a numerical experiment, it is way harder to have a precise value of cementation in Lab experiments, and the comparison with real fault gouge remains complicated. Based on the description of infill material (Riedmüller et al., 2001; Sibson, 1986; Wise et al., 1985), we can also consider another brecciation process for which the observed infill material is incohesive. This infill material is called “matrix” particles and refers to very thin particles produced by previous fragmentation or introduction of sediments within the fault gouge (Woodcock & Mort, 2008). The DEM approach seems also suitable for this kind of granular-type filling and consists of modelling a bi-disperse mixture between angular and faceted grains surrounded by hexagonal cells representing a matrix of fines. This represents a promising research avenue. Another improvement would consist in considering 3D DEM modelling in order to overcome the intrinsic limitations of 2D simulations, which can only claim to be very simplified conceptualization of the reality.

In future work, we also wish to extend this numerical campaign to cover a wider range of normal stresses and gouge thickness, in order to derive an empirical slip-weakening friction law based on the micromechanical properties of the gouge. This law could then be implemented in dynamic rupture modelling at a larger scale for a dialogue with seismological data. An interesting line of work could also be to investigate the time scale of the cementation evolution in a fault and to evaluate to what extent the associated evolution of the fault strength interacts with its seismic cycle.

### Acknowledgements

We thank Arnold Blaisonneau, Julie Maury, and Theophile Guillon from BRGM (France) for fruitful discussions all along with this research. We also thank François Passelègue for his kind remarks and discussions at the end of the paper. We express gratitude to the reviewing team for the interesting and instructive comments on the first version of the manuscript. This study was funded by INSA Lyon, by
means of a national PhD grant led by A. Daouadji and
G. Mollon.

Declarations

**Conflict of interest** The authors declare that they have no competing financial interests. Numerical data and software used for this research are available in these in-text data citation references: Casas et al. (2020), “Cohesion and Initial Porosity of Granular Fault Gouges”, Mendeley Data, V2, https://doi.org/10.17632/7c3dcj7spw.2. Other explanations are included in its supporting information file or available by contacting the corresponding author at nath27casas@gmail.com (Nathalie Casas).

REFERENCES

An, L.-J., & Sammis, C. (1994). Particle size distribution of cataclastic fault materials from southern California: a 3D study. *Pure and Applied Geophysics*, *143*(1), 203–227. https://doi.org/10.1007/BF00874329

Anthony, J. L., & Marone, C. (2005). Influence of particle characteristics on granular friction. *Journal of Geophysical Research: Solid Earth*, *110*(8), 1–14. https://doi.org/10.1029/2004JB003399

Antonellini, M. A., & Pollard, D. D. (1995). Distinct element modeling of deformation bands in sandstone. *Journal of Structural Geology*, *17*(8), 1165–1182. https://doi.org/10.1016/0191-8141(95)00001-T

Azéma, E., & Radjai, F. (2014). Internal structure of inertial granular flows. *Physical Review Letters*, *112*(7), 1–5. https://doi.org/10.1103/PhysRevLett.112.078001

Berger, N., Azéma, E., Douce, J. F., & Radjai, F. (2015). Scaling behaviour of cohesive granular flows. *EPL*. https://doi.org/10.1209/0295-5075/112/64004

Biegel, R. L., Sammis, C. G., & Dieterich, J. H. (1989). The frictional properties of a simulated gouge having a fractal particle distribution. *Journal of Structural Geology*, *11*(7), 827–846. https://doi.org/10.1016/0191-8141(89)90101-6

Billi, A. (2005). Grain size distribution and thickness of breccia and gouge zones from thin (< 1 m) strike-slip fault cores in limestone. *Journal of Structural Geology*, *27*(10), 1823–1837. https://doi.org/10.1016/j.jsgc.2005.05.013

Billi, A., & Storti, F. (2004). Fractal distribution of particle size in carbonate cataclastic rocks from the core of a regional strike-slip fault zone. *Tectonophysics*, *384*(1–4), 115–128. https://doi.org/10.1016/j.tecto.2004.03.015

Bishop, A. W. (1971). The influence of progressive failure on the choice of the method of stability analysis. *Geotechnique*, *21*(2), 168–172. https://doi.org/10.1680/geot.1971.21.2.168

Blenkinsop, T. G. (1991). Cataclasis and processes of particle size reduction. *Pure and Applied Geophysics*, *136*(1), 59–86. https://doi.org/10.1007/BF00878888

Byerlee, J. D., & Brace, W. F. (1968). Stick slip, stable sliding, and earthquakes—effect of rock type, pressure, strain rate, and stiffness. *Journal of Geophysical Research*, *73*(18), 6031–6037. https://doi.org/10.1029/JB073i018p06031

Byerlee, J. D., & Savage, J. C. (1992). Coulomb plasticity within the fault zone. *Geophysical Research Letters*, *19*(23), 2341–2344. https://doi.org/10.1029/92GL02370

Cho, N., Martin, C. D., & Sego, D. C. (2007). A clumped particle model for rock. *International Journal of Rock Mechanics and Mining Sciences*, *44*(7), 997–1010. https://doi.org/10.1016/j.ijrmms.2007.02.002

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

Cundall, P. A. (1989). Numerical experiments on localization in frictional materials. *Ingenieur-Archiv*, *59*, 148–159.

Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. *Geotechnique*, *29*(1), 47–65. https://doi.org/10.1680/geot.1979.29.1.47

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

Daouadji, A., & Hicher, P. (2010). An enhanced constitutive model for crushable granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics*, *34*, 555–580. https://doi.org/10.1002/nag

Das, A., Tengattini, A., Nguyen, G. D., Viggiani, G., Hall, S. A., & Einav, I. (2014). A thermomechanical constitutive model for cemented granular materials with quantifiable internal variables. *Part II—validation and localization analysis. Journal of the Mechanics and Physics of Solids*, *70*(1), 382–405. https://doi.org/10.1016/j.jmps.2014.05.022

Di Toro, G., Han, R., Hirose, T., Paola, N. D., Nielsen, S., Mizoguchi, K., et al. (2011). Fault lubrication during earthquakes. *Nature*, *471*, 494–498. https://doi.org/10.1038/nature09838

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

Dorostkar, O., Guyer, R. A., Johnson, P. A., Marone, C., & Carmeliet, J. (2017). On the role of fluids in stick-slip dynamics of saturated granular fault gouge using a coupled computational fluid dynamics-discrete element approach. *Journal of Geophysical Research: Solid Earth*, *122*(5), 3689–3700. https://doi.org/10.1002/2017JB014099

Dorostkar, O., Guyer, R. A., Johnson, P. A., Marone, C., & Carmeliet, J. (2018). Cohesion-induced stabilization in stick-slip dynamics of weakly wet, sheared granular fault gouge. *Journal of Geophysical Research: Solid Earth*, *123*(3), 2115–2126. https://doi.org/10.1002/2017JB015171

Dunn, D. E., LaFountain, L. J., & Jackson, R. E. (1973). Porosity dependence and mechanism of brittle fracture in sandstones. *Journal of Geophysical Research*. https://doi.org/10.1029/JB078i014p02403

Estrada, N., Lizcano, A., & Taboada, A. (2010). Simulation of cemented granular materials. I. Macroscopic stress–strain response and strain localization. *Physical Review E*, *82*, 1–11. https://doi.org/10.1103/PhysRevE.82.011303
Faqih, A., Chaudhuri, B., Muzzio, F. J., Tomassone, M. S., & Hammond, S. (2006). Flow-induced dilation of cohesive granular materials. *AIChE Journal, 52*(12), 4124-4132.

Ferdowsi, B. (2014). Discrete element modeling of triggered slip in faults with granular gouge application to dynamic earthquake triggering. *ETH Zurich*. https://doi.org/10.3929/ethz-a-010232908

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

Friedman, M., Handin, J., & Alani, G. (1972). Fracure-surface energy of rocks. *International Journal of Rock Mechanics and Mining Sciences and Geomechanics, 9*(6), 757-764. https://doi.org/10.1016/0148-9062(72)90034-4

Frye, K. M., & Marone, C. (2003). The effect of particle dimensionality on granular friction in laboratory shear zones. *Geophysical Research Letters, 29*(19), 22-1-22-4. https://doi.org/10.1029/2002gl015709

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

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

Handin, J. (1969). On the Coulomb-Mohr failure criterion. *Journal of Geophysical Research, 74*(22), 5343-5348. https://doi.org/10.1029/JB074i022p05343

Harmon, J. M., Karapiperis, K., Li, L., Moreland, S., & Andrade, J. E. (2021). Modeling connected granular media: Particle bonding within the level set discrete element method. *Computer Methods in Applied Mechanics and Engineering, 373*, 113486. https://doi.org/10.1016/j.cma.2020.113486

Hirth, G., & Tullis, J. (1989). The effects of pressure and porosity on the micromechanics of the brittle-ductile transition in quartzite. *Journal of Geophysical Research*. https://doi.org/10.1029/jb094ib12p17825

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

Jaeger, J. C. (1971). Friction of rocks and stability of slopes. *Geotechniques, 116*(4-5), 615-626. https://doi.org/10.1007/BF00876528

Jiang, M., Zhang, W., Sun, Y., & Utili, S. (2013). An investigation on loose cemented granular materials via DEM analyses. *Granular Matter, 15*(1), 65-84. https://doi.org/10.1007/s10035-012-0382-8

Kasyap, S. S., & Senetakis, K. (2020). An experimental investigation on the tribological behaviour of nominally flat quartz grains with gouge material in dry, partial saturated and submerged conditions. *Pure and Applied Geophysics*. https://doi.org/10.1007/s00024-020-02431-1

Kasyap, S., & Senetakis, K. (2021). Micromechanical-based experimental and analytical studies on rate effects and stick-slip instability of smooth quartz surfaces in the presence of plastic and non-plastic gouges. *International Journal for Numerical and Analytical Methods in Geomechanics, 45*(6), 776-793. https://doi.org/10.1002/nag.3178

Kranz, R. L., & Scholz, C. H. (1977). Critical dilatant volume of rocks at the onset of tertiary creep crack observation using the scanning electron microscope each quartzite. *Journal of Geophysical Research*. https://doi.org/10.1029/JB082i030p04893

Lade, P. V., & Overton, D. D. (1989). Cementation effects in frictional materials. *Journal of Geotechnic Engineering, 115*(10), 1373-1387. https://doi.org/10.1061/(ASCE)0733-9410(1989)115:10(1373)

Leeman, J., Scuderi, M. M., Marone, C., & Saffer, D. (2015). Stiffness evolution of granular layers and the origin of repetitive, slow, stick-slip frictional sliding. *Granular Matter, 17*(4), 447-457. https://doi.org/10.1007/s10035-015-0565-1

Lin, A. (1999). S-C cataclasite in granitic rock. *Tectonophysics, 304*(3), 257-273. https://doi.org/10.1016/S0040-1951(99)00026-8

Lockner, D. A., & Beeler, N. M. (2002). Rock failure and earthquakes. *International Geophysics Series, 81*(A), 505-538.

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

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

Majmudar, T. S., & Behringer, R. P. (2005). Contact force measurements and stress-induced anisotropy in granular materials. *Nature, 435*(June), 1079-1082. https://doi.org/10.1038/nature03805

Marone, C. (1995). Fault zone strength and failure criteria. *Geophysical Research Letters, 22*(6), 723-726. https://doi.org/10.1029/95GL00268

Marone, C. (1998). Laboratory-derived friction laws and their application to seismic faulting. *Annual Review of Earth Planetary Sciences, 26*, 643-696.

Marone, C., Hobbs, B. E., & Ord, A. (1992). Coulomb constitutive laws for friction: Contrasts in frictional behavior for distributed and localized shear. *Pageoph. https://doi.org/10.1007/BF00876327*

Marone, C., & Kilgore, B. (1993). Scaling of the critical slip distance for seismic faulting with shear strain in fault zones. *Letters to Nature, 362*, 210-211. https://doi.org/10.1038/362618a0

Marone, C., & Scholz, C. H. (1988). The depth of seismic faulting and the upper transition from stable to unstable slip regimes. *Geophysical Research Letters, 15*(8), 621-624. https://doi.org/10.1029/GL015i006p00621

Marone, C., & Scholz, C. H. (1989). Particle-size distribution and microstructures within simulated fault gouge. *Journal of Structural Geology, 11*(7), 799-814. https://doi.org/10.1016/0191-8141(89)90099-0

Maurer, W. C. (1965). Shear failure of rock under compression. *SPE Journal, 5*, 167-176. https://doi.org/10.2118/1054-PA

Mead, W. J. (1925). The geologic role of dilatancy. *The Journal of Geology, 33*(7), 685-698. https://doi.org/10.1086/623241

Meng, F., Wong, L. N. Y., & Zhou, H. (2021). Rock brittleness indices and their applications to different fields of rock engineering: A review. *Journal of Rock Mechanics and Geotechnical Engineering, 13*(1), 221-247. https://doi.org/10.1016/j.jrmge.2020.06.008
Midi, G. R. D. (2004). On dense granular flows. *European Physical Journal E, 14*(4), 341–365. https://doi.org/10.1140/epje/i2003-10153-0

Misra, S., Mandal, N., & Chakraborty, C. (2009). Formation of Riedel shear fractures in granular materials: Findings from analogue shear experiments and theoretical analyses. *Tectonophysics, 471*(3–4), 253–259. https://doi.org/10.1016/j.tecto.2009.02.017

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. (2018a). A unified numerical framework for rigid and compliant granular materials. *Computational Particle Mechanics, 5*, 517–527. https://doi.org/10.1007/s40571-018-0187-6

Mollon, G. (2018b). Mixtures of hard and soft grains: Micromechanical behavior at large strains. *Granular Matter*. https://doi.org/10.1007/s10035-018-0812-3

Mollon, G., Quacquarelli, A., Andò, E., & Viggiani, G. (2020). Can friction replace roughness in the numerical simulation of granular materials? *Granular Matter*. https://doi.org/10.1007/s10035-020-1004-5

Mollon, G., & Zhao, J. (2012). Fourier–Voronoi-based generation of realistic samples for discrete modelling of granular materials. *Granular Matter, 5*, 621–638. https://doi.org/10.1007/s10035-012-0356-x

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

Morgan, J. K., & Boettcher, M. S. (1999). Numerical simulations of granular shear zones using the distinct element method: 1. Shear zone kinematics and the micromechanics of localization. *Journal of Geophysical Research: Solid Earth, 104*(B2), 2703–2719. https://doi.org/10.1029/1998jb900056

Muto, J., Nakatani, T., Nishikawa, O., & Nagahama, H. (2015). Fractal particle size distribution of pulverized fault rocks as a function of distance from the fault core. *Geophysical Research Letters, 42*(10), 3811–3819. https://doi.org/10.1002/2015GL064026

Niemeijer, A., Toro, G. D., Nielsen, S., & Felice, F. D. (2011). Frictional melting of gabbro under extreme experimental conditions of normal stress, acceleration, and sliding velocity. *Journal of Geophysical Research*. https://doi.org/10.1029/2010JB008181

Noel, C., Passelègue, F. X., & Violay, M. (2021). Brittle faulting of ductile rock induced by pore fluid pressure build-up. *Journal of Geophysical Research*. https://doi.org/10.1029/2020JB021331

Nouguier-lehon, C., Cambou, B., & Vincens, E. (2003). Influence of particle shape and angularity on the behaviour of granular materials: A numerical analysis. *International Journal for Numerical and Analytical Methods in Geomechanics, 27*(14), 1207–1226. https://doi.org/10.1002/nag.314

Olgaard, D. L., & Brace, W. F. (1983). The microstructure of gouge from a mining-induced seismic shear zone. *International Journal of Rock Mechanics and Mining Sciences and Geomechanics, 20*(1), 11–19. https://doi.org/10.1016/0148-9062(83)91610-8

Paglialunga, F., Passelègue, F., Brantut, N., Barras, F., Lebihain, M., & Violay, M. (2021). On the scale dependence in the

dynamics of frictional rupture: Constant fracture energy versus size-dependent breakdown work. *Preprint in Earth and Planetary Science Letters*, 1–33. arXiv:2104.15103

Philit, S., Soliva, R., Castilla, R., Ballas, G., & Taillefer, A. (2018). Clusters of cataclastic deformation bands in porous sandstones. *Journal of Structural Geology, 114*, 235–250. https://doi.org/10.1016/j.jsg.2018.04.013

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

Rathbun, A. P., Renard, F., & Abe, S. (2013). Numerical investigation of the interplay between wall geometry and friction in granular fault gouge. *Journal of Geophysical Research: Solid Earth, 118*(3), 878–896. https://doi.org/10.1002/jgrb.50106

Rice, J. (2006). Heating and weakening of faults during earthquake slip. *Journal of Geophysical Research: Solid Earth, 111*(5), 1–29. https://doi.org/10.1029/2005JB004006

Rice, J., & Cocco, M. (2002). Seismic fault rheology and earthquake dynamics. *Journal of Geophysical Research, 107*, 1000–1029.

Riedmüller, B. G., Brosch, F. J., Klima, K., & Medley, E. W. (2001). Engineering geological characterization of brittle faults and classification of fault rocks. *Engineering Geology, 19*(April), 13–19.

Rognon, P. G., Roux, J. N., Naaim, M., & Chevoir, F. (2008). Dense flows of cohesive granular materials. *Journal of Fluid Mechanics, 596*, 21–47. https://doi.org/10.1017/S0022112007009329

Roy, S., & Luding, S. (2017). Effect of cohesion on local compaction and granulation of sheared soft granular materials. *Powders and Grains, 03065*, 2–5. https://doi.org/10.1051/epjconf/201714003065

Sammis, C., & Biegel, R. (1989). Fractals, fault-gouge, and friction. *Pure and Applied Geophysics, 131*, 255–271. https://doi.org/10.1007/BF00874490

Sammis, C., King, G., & Biegel, R. (1987). The kinematics of gouge deformation. *Pure and Applied Geophysics, 125*, 777–812. https://doi.org/10.1007/BF00878033

Sandeep, C. S., & Senetakis, K. (2019). An experimental investigation of the microslip displacement of geological materials. *Computers and Geotechnics, 107*(June 2018), 55–67. https://doi.org/10.1016/j.compgeo.2018.11.013

Schellart, W. P. (2000). Shear test results for cohesion and friction coefficients for different granular materials: Scaling implications for their usage in analogue modelling. *Tectonophysics, 324*(1–2), 1–16. https://doi.org/10.1016/S0040-1951(00)00111-6

Sibson, R. H. (1986). Brecciation processes in fault zones: Inferences from earthquake rupturing. *Pure and Applied Geophysics PAGEOPH, 124*(1–2), 159–175. https://doi.org/10.1007/BF00875724

Taylor, D. W. (1948). *Fundamentals of soil mechanics*. (J. W. & Sons, Ed.). London: CHAPMAN & HALL, Limited.

Tchalenko, J. S. (1970). Similarities between shear zones of different magnitudes. *Geological Society of America Bulletin, 81*(6), 1625–1640. https://doi.org/10.1130/0016-7606(1970)81[1625:SBSZOD]2.0.CO;2

Tengattini, A., Das, A., Nguyen, G. D., Viggiani, G., Hall, S. A., & Einav, I. (2014). A thermomechanical constitutive model for cemented granular materials with quantifiable internal variables.
Part I—theory. *Journal of the Mechanics and Physics of Solids*, **70**(1), 281–296. https://doi.org/10.1016/j.jmps.2014.05.021

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

Wise, D. U., Dunn, D. E., Engelder, T., Geiser, P. A., Hatcher, R. D., Kish, S. A., et al. (1985). Fault-related rocks: Suggestions for terminology. *Geology*, **13**(3), 218. https://doi.org/10.1130/0091-7613(1985)13%3c218b:CAROFB%3e2.0.CO;2

Wissa, A. E. Z. (1965). *Effective stress-strength behavior of cemented soils (Thesis (Sc D.)).* Massachusetts Institute of Technology.

Wood, D. M. (1990). Soil behaviour and critical state soil mechanics. *Cambridge University Press, Cambridge*. https://doi.org/10.1201/9781351255400-1

Woodcock, N. H., & Mort, K. (2008). Classification of fault breccias and related fault rocks. *Geological Magazine*, **145**(3), 435–440. https://doi.org/10.1017/S0016756808004883

Zhang, L., Nguyen, N. G. H., Lambert, S., Nicot, F., Prunier, F., & Djeran-Maigre, I. (2017). The role of force chains in granular materials: From statics to dynamics. *European Journal of Environmental and Civil Engineering*, **21**(7–8), 874–895. https://doi.org/10.1080/19648189.2016.1194332

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