# Self-affine roughness and spatial-temporal correlation of shear-induced micro-fractures in rough DEM rock joints

Alberto Varela-Valdez, Gilles Pijaudier-Cabot, Moisés Hinojosa Rivera, Victor

I German-Flores, Stéphane Morel

► To cite this version:

Alberto Varela-Valdez, Gilles Pijaudier-Cabot, Moisés Hinojosa Rivera, Victor I German-Flores, Stéphane Morel. Self-affine roughness and spatial-temporal correlation of shear-induced micro-fractures in rough DEM rock joints. Engineering Fracture Mechanics, 2024, pp.110767. 10.1016/j.engfracmech.2024.110767. hal-04859176

HAL Id: hal-04859176

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

Submitted on 30 Dec 2024

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

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

<figure><img src="image_1.png" /></figure>
# Self-Affine Roughness and Spatial-Temporal Correlation of Shear-Induced Micro-Fractures in Rough DEM Rock Joints

Alberto Varela Valdez<fnref n="1" />, Gilles Pijaudier-Cabot<fnref n="2" />, Moisés Hinojosa Rivera<fnref n="3" />, Victor I. German-Flores<fnref n="3" />, Stéphane Morel<fnref n="4" />.

## Abstract

We investigate the spatial and temporal correlations of micro-fractures generated during shear simulations of rough rock joints using the Discrete Element Method (DEM). Eight numerical joints were meticulously designed, characterized by distinct self-affine roughness parameters: the Hurst exponent *H*, self-affine correlation length *L*<sub>c</sub>, and variance of heights *σ*<sup>2</sup>. Shearing was conducted under quasi-static conditions, maintaining constant normal load (CNL) and constant shear speed. The DEM model underwent calibration to faithfully replicate the mechanical behavior observed in experimental compression, tension and shear tests on synthetic rock. Remarkably, our results unveil a power-law behavior in both spatial and temporal dimensions. The probability distribution *P*(*Δd*) describing the separation distance *Δd* between two fractures follows a power-law relationship: *P*(*Δd*)~*Δd*<sup>-q</sup>, where the exponent *q* ranges between 0.4 and 0.8, displaying a marked sensitivity to joint roughness. Additionally, the probability *p*(*Δt*) of two fractures occurring within a time interval *Δt* during shearing also conforms to a power-law distribution: *P*(*Δt*)~*Δt*<sup>-p</sup>. For mechanically interacting micro-fractures, *n* ranges from 0.716 to 0.869, while for micro-fractures generated across the entire joint area (not necessarily interacting fractures), *n* takes values between 1.458 and 1.895. Notably, these findings underscore a remarkable analogy with the Omori law in seismic activity.

**Keywords:** Rock joints, Discrete Element Method, Self-Affine Roughness, Direct Shear, Friction, Micro-cracks

## 1. Introduction

It is well known in rock mechanics and reservoir geo-mechanics, that pre-existing fractures (or joints) plays a major role on the thermo-hydro-mechanical behavior of the rock masses [¹⁻⁵]. One of the key parameters that controls such behavior is the roughness of rock joints [⁴,⁶⁻⁹], which exhibits a self-affine regime that expands on several order of magnitudes [¹⁰,¹¹]. One characteristic of the self-affine profiles

<footnote n="1">¹Departamento de Minería, Instituto Tecnológico Superior Zona Occidente, Av. Tecnológico 2000, Loma la Perla 99102 Sombrerete, Zacatecas, México</footnote>

<footnote n="2">²Université de Pau et des Pays de l'Adour, CNRS, TotalEnergies, LFCR, Allée du parc Montaury, 64600 Anglet, France.</footnote>

<footnote n="3">³Facultad de Ingeniería Mecánica y Eléctrica, Universidad Autónoma de Nuevo León, Av. Universidad s/n, Ciudad Universitaria, 66451 San Nicolás de los Garza, Nuevo León, México.</footnote>

<footnote n="4">⁴Université de Bordeaux, UMR 5295, Institut de Mécanique et d'Ingénierie - Bordeaux (I2M), Département Génie Civil et Environnemental (GCE), Bordeaux, France.</footnote>
is that the mean height difference $\Delta h(\varepsilon) = \langle [h_{r+\varepsilon} - h_r]^2 \rangle$ for two points on the surface separated by a distance $\varepsilon \ll L_c$ follows a power law of the form $\Delta h(\varepsilon) = B\varepsilon^{2(1-H)}$, where H is the Hurst exponent, whose values range from 0.4 to 0.8 [12–14] and $L_c$ is the self-affine correlation length and B constant (called topothesy). For $\varepsilon \gg L_c$, the rough profile is no longer self-affine and $\Delta h(\varepsilon \gg L_c) = 2\sigma^2$, where $\sigma^2$ is the variance of heights of the rough surface.

The dependence of mechanical and hydro-mechanical behavior on the self-affine fracture roughness has been studied by several authors [15–18], in particular, special attention has been paid on the dependency between the Hurst exponent H, the shear stress τ and the dilatancy δ of the rock joint [19–21]. It has been reported that as the Hurst exponent decreases, the dilatancy and the shear stress peak $\tau_{pic}$ increase [19,21]. The variance of heights $\sigma^2$ of the rough fracture, as well the self-affine correlation length $L_c$, also influence the shear behavior of the joints [19] i.e., as $\sigma^2$ increases, $\tau_{pic}$ increases, and as $L_c$ increases, $\tau_{pic}$ decreases [19]. Also, it is common that channels are formed perpendicular to the shear direction, which results on an increase of the hydraulic permeability of the fracture [19,22,23].

During the sliding of rock joints, acoustic emission is produced because of the micro-fracturing of the rock [24,25]. A peak of the acoustic emission events is obtained as the shear stress approaches to $\tau_{pic}$ [26,27]. The statistics of acoustic emissions (AE) in intact rock has been studied by several authors, e.g. by Rudajev et al. [28]. They found that the AE of sandstone subjected to compression load follows a Omori's law type of behavior, $n(t) \sim t^{-p}$, where n(t) is the number of aftershocks occurred after the main quake, t is the time passed after the main quake and p is a parameter taking values between 0.5-3.0 [29]. Rudajev et al. observed also that p decreases as the compression stress increases. Bunger et al. [30] adjusted a modified Omori's formula to AE from hydraulic fracture test on granite, the best fit was done taking $p = 1$. Ojala et al. [31] carried out confined compression test at different temperatures and strain rates on porous sandstone. The AE emitted during the test follows a Omori's law with the p-values related linearly with temperature. The p-values reported by Ojala et al. varied from 0.79 to 1.97 for the aftershocks and from 0.80 to 6.81 for foreshocks. Smirnov et al. [32] carried out confined compression tests on intact and fractured specimens of sandstone. They found that p increases with compression stress for the fractured specimens, but not change on p was noticed for the intact specimens as the compression stress increases.

The above statistical analyses are related to intact rock specimens, not to rock joints. It remains unclear how the acoustic emission or micro-fracturing is related to the roughness during the sliding of the rock joints. Thus, the main goal of this work is to present numerical analyses that should help at understanding better the dependency between self-affine joint roughness and micro-fracturing during the shearing of rock joints.

This paper is organized as follows: section 2 provides a recap of the numerical discrete model utilized for the analysis and outlines the methodology for its calibration. In Section 3, temporal correlations in micro-fractures induced by joint shear are discussed, while Section 4 deals with spatial correlations. Finally, Section 5 is a summary of the primary findings and conclusions.
## 2. Methodology

The present analysis relies on a description of rock joints with the discrete element method (DEM) that has been already implemented in a previous study [¹⁹]. The joint is considered to possess a given thickness and the main steps of this modelling methodology are recalled in the following.

First, we start with the calibration of a representative statistical volume element (SEV) that depicts the mechanical response of the joint. For this purpose, we consider a 1 cm side cube filled with spherical Discrete Elements (DEs). The sizes of these DEs follow a Gaussian distribution with a mean of 0.5 mm and a standard deviation of 0.1 mm. The DEs, with a density of 2650 kg/m³, normal stiffness of 1 N/m, tangential stiffness of 0.1 N/m, and a friction coefficient of 0.5, adhered to Newton's laws of motion. Interconnections between DEs are established to form a solid. The connections are cylindrical beams with a diameter equal to the smallest DE being joined. These connections, referred to as cylindrical or parallel bonds according to the nomenclature used by PFC3D of Itasca (the software employed for DEM simulations), are considered continuous, linear-elastic, and follow an elliptical fracture criterion. When the maximum normal and shear stresses meet the rupture criterion, these bonds break, enabling the simulation of the macroscopic behavior of quasi-brittle materials like concrete, mortar, or rocks. Additionally, it's worth noting that the parallel bonds (cylindrical beams) also have normal and tangential stiffness of 0.1 N/m and 0.01 N/m, respectively. The detailed parameters for this calibration, including elasticity modulus, Poisson's coefficient, maximum compression strength, and maximum tensile strength, are outlined in Table 1.

<table><caption>Table 1. Comparison of DEM Model with experimental results and roughness parameters</caption><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>S1</td><td>H =0.53, Lc = 1.6 cm, σ² = 0.547 mm²</td><td rowspan="2">DEM</td><td rowspan="2">Experimental data</td></tr><tr><td>S2</td><td>H =0.53, Lc = 1.6 cm, σ² = 0.328 mm²</td></tr><tr><td>S3</td><td>H =0.53, Lc = 2.0 cm, σ² = 0.547 mm²</td><td rowspan="2">Compressive strength</td><td>71 MPa</td></tr><tr><td>S4</td><td>H =0.53, Lc = 2.0 cm, σ² = 0.328 mm²</td><td>74.5 MPa</td></tr><tr><td>S5</td><td>H =0.41, Lc = 1.6 cm, σ² = 0.547 mm²</td><td>Tension strength</td><td>7.0 MPa</td></tr><tr><td>S6</td><td>H =0.41, Lc = 1.6 cm, σ² = 0.328 mm²</td><td>Poisson Coefficient</td><td>0.19</td></tr><tr><td>S7</td><td>H =0.41, Lc = 2.0 cm, σ² = 0.547 mm²</td><td>Young Modulus</td><td>32 MPa</td></tr><tr><td>S8</td><td>H =0.41, Lc = 2.0 cm, σ² = 0.328 mm²</td><td></td><td>31 MPa</td></tr></tbody></table>

Second, once the material describing the joint has been properly described with DEM, we turn to the generation of the joints with a self-affine roughness. We consider a joint size of 10 cm in length and width, with a height of 1 cm (matching the size of the SEV). Approximately 100,000 Discrete Elements (DEs), possessing the same characteristics as those used in the calibration are used to construct these joints. The process of creating joints with controlled roughness involves the initial generation of numerical joint molds in which the DEs are placed. These molds, each exhibiting the desired roughness, are crafted using the power spectral method [¹⁹]. Through this method, self-affine rough surfaces can be generated while precisely controlling the variance of heights σ², the correlation length Lc, and the roughness exponent H. This process ensures the replication of desired roughness characteristics in the generated rock joints. The details of the roughness parameters manipulated during this process are outlined in Table 1. It is important to recall that all parameters were rigorously validated against experimental roughness data [¹⁹]. Figure 0 illustrates one of the obtained joint.
<figure><img src="image_2.png" /><figcaption>Figure 0 DEM joint</figcaption></figure>

Finally, we may advance to conduct friction tests under constant normal load conditions on the rock joint. A constant normal load in the direction of the z-axis, set at 14 MPa, is applied initially. Then, the upper part of the joint is subjected to a constant sliding velocity of 6x10⁻⁵ m/s, while the lower part remains fixed. Throughout these numerical tests, friction forces, vertical displacement of the joint, spatial localization of micro-fractures, and the timing of their occurrence are recorded.

Figure 1 presents the shear response and the evolution of the dilatancy upon shearing that was obtained in Ref. [¹⁹]. Throughout this paper, we will place particular emphasis on the joints DEMS4 and DEMS5 which are the smoothest and roughest joints respectively. However, the analysis covers all the joints considered in Table 1 whose response is depicted in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1. Shear response (left) and dilatancy (right) of the considered joints (after ref. 19).</figcaption></figure>

Having successfully obtained the mechanical response of our DEM model, we examine now the micro-fracturing process, more precisely temporal and spatial correlations, with the aim to elucidate relationships between this process and the self-affine roughness characteristics.

### 3. Temporal correlation of micro-fractures
This section investigates the possible temporal correlations between micro-fractures that were created during the shearing of DEM joints. It is useful to point out that shearing induces a complex fracture process: when a bond breaks (micro-fracture), part of the elastic energy stored in the bonds is dissipated in the creation of a micro-fracture. The other part is converted into kinetic energy that may trigger subsequent micro-fractures as part of a stress redistribution process. On one hand, once a micro-fracture is created, more fractures around the first one [³³] in a small-time interval Δt, can be formed. This kind of fractures are going to be named as interacting micro-fractures. On the other hand, during shearing, the two surfaces of the joint can be in contact on several points apart due to the roughness. If the stresses at these contact points are intense, micro-fractures may occur at locations separated by large distances. This kind of fractures are denoted as non-interacting micro-fractures.

Estimating the distance or time interval from which the interactions between an existing micro-fracture and subsequent ones do not interact is a challenging issue that requires to distinguish crack interactions from load redistribution. Theoretically, in a full dynamic process, redistribution occurs at a speed that is in relationship with the wave velocity. Within a sphere centered at the location of a given micro-fracture denoted as F¹, subsequent micro-fractures F¹ that occur at a time interval that is greater than the time needed for the wave to propagate from F¹ to F¹ are assumed to result from redistribution. It is considered that they are not interacting, whereas the remaining micro-fractures observed in the sphere are interacting ones.

For the sake of discussion, let us assume that the mean speed of propagation, $v_p$ of the elastic P-waves in the DEM model can be approximated as. [³⁴]:

$$v_p \approx D \sqrt{\frac{k(1-\alpha)}{m}} \qquad (1)$$

where D is the mean diameter of the particles, k is the stiffness of the bonds, α is the local damping coefficient and m is the mean mass of the particles. Introducing the corresponding values in Eq. (1), a mean value for the P-wave velocity $v_p \approx 9.8 \text{ mm/s}$ is obtained. Then the time interval limit between interacting and non-interacting subsequent micro-fractures will be set arbitrarily. It will correspond to the time needed for the wave to propagate through the height of the joint.

### 3.1. Interacting fractures produced during shearing

Figure 2*Érreur ! Source du renvoi introuvable.a* shows all the fractures (black dots) that were produced after 6mm of shearing of the DEMS5 joint. Consider a small set of fractures F¹ (i=1: N, N = number of fractures) formed during a time period *to* + Δ*t*<sub>p</sub> and localized inside a sphere of radius Δ*d*<sub>p</sub> centered at *do* (Figure 2b), *to* is the time at which the first fracture F¹ is formed and *do* is its position. The value of Δ*d*<sub>p</sub> is fixed at 10 mm, the height of the DEM joint, and Δ*t*<sub>p</sub> is the time that the elastic P-waves takes to travel the distance Δ*d*<sub>p</sub> (10 mm) at a velocity *v*<sub>p</sub>.
<figure><img src="image_4.png" /></figure>

**Figure 2** Schematic representations for the calculation of waiting times between fractures represented by small circles.

Then, the differences in time $\Delta t$ (or waiting times) between the formation of fracture $F^1$ and the rest of fractures $F^i$ occurring inside the sphere of radius $\Delta d_p$ is computed. Because all the fractures $F^i$ are created at time intervals $\Delta t$ smaller that the time needed of the P-wave to travel the distance $\Delta d_p$, these are considered as interacting fractures.

As an example, Figure 3 shows the probability distribution, i.e., histogram formed from the waiting times $\Delta t$, for the joints DEM S4 and DEM S5. This figure reflects the possible interactions between fractures that are localized in time and space.

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

**Figure 3** Distribution of waiting times $\Delta t$ between fractures. $p=0.869$ for DEMS4 (blue curve) and $p=0.720$ for DEMS5 (red curve).

In this figure, $P_{2f}(\Delta t)$ can be interpreted as the conditional probability that a fracture $F^i$ is formed after a time interval $\Delta t$ once the fracture $F^1$ is created. It can be seen from Figure 3 that for small time intervals ($\Delta t < 10^{-2}$ s) there is relatively high probability that a secondary fracture occurs, then, the
conditional probability for the formation of the second fracture decreases following a power law behavior reminiscent of Omori's law for the aftershocks [²⁹]:

$$P_{2f}(\Delta t) = \frac{A}{(\Delta t)^p} \qquad (2)$$

where A, and p are constants. The exponent p estimated from the fitted equation to the numerical data in Figure 3 ranges between 0.716-0.869 (Table 2), similar to the experimental values for the foreshocks and aftershocks reported for earthquakes [²⁹,³⁵], for the fracture of sandstone in confined compression test [³¹,³²], and for earthquakes produced by volcanic eruptions [³⁶]. Note that the exponent p depends slightly on the joint roughness: when the roughness of the joint decreases (DEM S4) p increases and vice versa for the roughest joint DEM S5. Finally, after a waiting time of around 0.11 seconds, the conditional probability of occurrence of the two fractures takes a relatively constant low value (Figure 3), indicating that the fracture process is not correlated in time afterwards.

## 3.2. Interacting and noninteracting fractures produced during shearing

As shear proceeds, successive fractures are formed (broken bonds), where each fracture $F^i$, is formed at a time $t^i$. To study the possible time-correlation between consecutive fractures $F^i$ and $F^{i+1}$, whether they interact or not, the differences in time $\Delta t = t^{i+1} - t^i$ between the formation of both fractures is computed. Figure 4 shows the probability distribution of the waiting times $\Delta t$ for all the fractures as they were produced chronologically during the shear simulations. The waiting time distribution for all fractures also follows a power law behavior given by Eq. (2), but with exponents p that ranges between 1.458 to 1.895, which seems to be systematically higher than the p exponent obtained for interacting fractures (p = 0.716-0.869).

<figure><img src="image_6.png" /><figcaption>Figure 4 Waiting time distribution for all the fractures and for joints DEMS4 and DEMS5. $p_{DEMS5} = 1.895$ (red curve), $p_{DEMS4} = 1.458$ (blue curve).</figcaption></figure>
Recall that carrying out such an analysis implies to consider the waiting times $\Delta t$ between fractures that not necessarily lies in the region enclosed by the red circle of **错! Source du renvoi introuvable.**, but also far away from each other. This is the reason why the p exponents in Figure 4 are different from those in Figure 3. Table 2 summarizes the results for all the 8 DEM joints.

It is important to remark that in Figure 4, the steady-state regime (as shown in Figure 3 with constant probability) is not observed. This is because the shear displacement is only 6 mm, whereas the minimum correlation length of the joints is 16 mm. There is a strong correlation between the upper and lower parts of the joint through the roughness. This will be observed in Figure 7, where micro-fractures concentrate in specific zones. For uncorrelated frictional surfaces, homogeneous micro-fracturing would be expected across the surface.

It is intriguing to observe the power law $P(t) \sim t^{-p}$ has been observed not only in foreshocks and aftershocks [29] but also in Acoustic Emissions (AE) during confined compression tests [31,32] and in the fracture of 2D materials like paper subjected to tension [37], suggesting a potential universality in fracture behavior.

<table><caption>Table 2 Exponents of the waiting time distribution between fractures</caption><thead><tr><th>Joint</th><th>p<sub>1</sub> (interacting fractures)</th><th>P<sub>2</sub> (all fractures)</th></tr></thead><tbody><tr><td>S1</td><td>0.810</td><td>1.616</td></tr><tr><td>S2</td><td>0.802</td><td>1.480</td></tr><tr><td>S3</td><td>0.831</td><td>1.481</td></tr><tr><td>S4</td><td>0.869</td><td>1.458</td></tr><tr><td>S5</td><td>0.720</td><td>1.895</td></tr><tr><td>S6</td><td>0.785</td><td>1.653</td></tr><tr><td>S7</td><td>0.716</td><td>1.557</td></tr><tr><td>S8</td><td>0.826</td><td>1.690</td></tr><tr><td>average</td><td>0.794</td><td>1.604</td></tr></tbody></table>

As explained earlier, for a self-affine profile, the structure of the surface follows the relationship, $\Delta h(\epsilon) = B\epsilon^{2(1-H)}$, when $\epsilon = L_c$, $\Delta h(\epsilon) = 2\sigma^2$, so the topothecy B of the self-affine surface can be expressed as $B = 2\sigma^2/L_c^{2(1-H)}$. Figure 5 illustrates the dependence of the exponents $p_1$ and $p_2$ observed for interacting micro-fractures and for all fractures respectively, on the surface roughness expressed in terms of the topothecy B.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Dependence of the p-exponent with the roughness of the joint expressed in terms of tis topothesis B. Exponents $p_1$ and $p_2$ correspond to interacting micro-fractures and to all fractures respectively.</figcaption></figure>

On this figure, it can be observed that the roughness of the joint has a significant effect on the micro-fractures generated during friction. There is a linear relationship between the roughness and the exponent $p_2$ for all fractures with slope $b_2 \approx 25$ and there is another linear relationship in the case of interacting fractures with slope $b_1 \approx -8.7$. Fractures linked to p2 are considered to do not interact through elastic waves but through the joint roughness. In fact, exponent $p_1$ might be considered as being almost independent from the roughness of the joint. The range of variation is small compared to that of exponent $p_2$. Micro-fractures that occur close to each other (interacting micro-fracture linked to p1) are not sensitive to the roughness which is more influent when micro-fractures are far apart (fractures linked to p2).

**4. Spatial correlation between fractures**

Figure 6 and Figure 7 show the spatial distribution of the fractures at the end of the shear simulation in the joints DEM S5 and DEM S4 (the joints were sheared in the –X direction).
<figure><img src="image_8.png" /><figcaption>Figure 6. Distribution of fractures at the end of the shear simulations for the joint DEM S5.</figcaption></figure>

The fractures form clusters in the X-Y plane that normally corresponds to the location of the higher asperities in the joints, see for instance (Figure 8) [19]. Further, the fractures tend to be aligned in bands in the planes X-Z and Y-Z (zones were the density of fractures increases and zones where there are almost not fractures).

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7. Distribution of fractures at the end of the shear simulations for the joint DEM S4.</figcaption></figure>

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

<figure><img src="image_12.png" /><figcaption>Figure 8. Localization of micro fractures (red dots) after shearing 6mm of the joint DEM S4</figcaption></figure>

For the sake of qualitative comparison, Figure 9 shows the distribution of earthquakes along an active fault in Turkey as reported in [³⁸]. Note that at the planetary level seems that the seismic zones also tend to be localized in bands perpendicular to the fault. Figure 9 seems to be qualitatively like the present results on small joint (although there is, without doubt some effects of the scale and accuracy of the measurement of the seismicity).
<figure><img src="image_13.png" /></figure>

Figure 9. Seismicity in the Çınarcık basin since 17 August 1999 (yellow: 17 August 1999–12 November 1999, red: 2001–2010) from [³⁸].

Figure 10 shows the probability $P_{2f}$ to find two fractures separated by a distance $\Delta d$. The probability to find two fractures at small distances is high and constant. Then, it decreases following a power law $P_{2f}(\Delta d) \sim \Delta d^{-q}$.

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

Figure 10. Distribution of distance $\Delta d$ between fractures, $q=0.84$ for DEMS4 (blue) and $q=0.46$ for DEMS5 (red).

The $q$ exponents obtained from linear regression range from 0.4 to nearly 0.85, with joint DEM S6 exhibiting the lowest value of 0.4, and joint DEM S1 displaying the highest value. The power-law relationship $P_{2f}(\Delta d) \propto \Delta d^{-q}$ can be integrated to obtain the spatial correlation function [³⁹], denoted as $C(\Delta d, w)$, which represents the normalized sum of fracture pairs with distances less than $\Delta d$ and $w$ is the moment of $C$.
If the set of distances between these fracture pairs forms a fractal object, then $C(\Delta d, w)$ should be proportional to $\Delta d^D$, where $D = -q + 1$ represents the fractal dimension of the set off distances. Indeed, considering a fractal object with dimension $D = 1 - H$, where $H$ represents the Hurst exponent of the fractal set of fractures, it follows that the $q$ exponent is equivalent to the Hurst exponent $H$. Therefore, by averaging the $q$ values (see Table 3), the fractal dimension $D$ of the distance set between fractures is approximately $D \approx 0.4$. Similar $q$ ($H$ values) values of 0.4, and 0.8 have also been observed in earthquake epicenter distances [⁴⁰].

<table><caption>Table 3. Localized fractures in time and space (coupled and uncoupled fractures)</caption><thead><tr><th></th><th>q</th></tr></thead><tbody><tr><td>S1</td><td>0.894</td></tr><tr><td>S2</td><td>0.624</td></tr><tr><td>S3</td><td>0.70</td></tr><tr><td>S4</td><td>0.843</td></tr><tr><td>S5</td><td>0.459</td></tr><tr><td>S6</td><td>0.400</td></tr><tr><td>S7</td><td>0.451</td></tr><tr><td>S8</td><td>0.409</td></tr><tr><td>average</td><td>0.60</td></tr></tbody></table>

The process of friction gives rise to distinct regions of fracture clouds that span into the joint from the frictional surface along the $\pm z$ direction (see Figs. 6,7). Within each of these observed fracture clouds, it is conceivable that a macroscopic fracture is forming. In other words, these micro-fracture clusters can be regarded as a process zone, as experimentally detected through acoustic emissions [⁴¹]. Ponson et al. [⁴²] studied the self-affine behavior of fractures induced by indirect tension in sandstone. They showed that fracture surfaces exhibit self-affine behavior characterized by a roughness exponent, $\zeta$, ranging from 0.4 to 0.5. Similarly, Morel et al. [⁴³] found similar roughness exponents, $\zeta = 0.4$ and $\zeta = 0.8$, for fracture surfaces in mortar. The self-affine nature of these surfaces implies that the height difference, $\delta z(dx)$, between two points separated by a distance $dx$ in the rough surfaces follows a power-law $\delta z(dx) \sim dx^\zeta$.

As discussed above, the micro-fracture cloud exhibits fractal characteristics, and within this fracture cloud, macroscopic cracks form and propagate. These macro-cracks also demonstrate a self-affine behavior [⁴²,⁴³]. It might then be proposed that the self-affine behavior of these macro-cracks emerges from their formation and propagation within a fractal micro-fracture cloud. In this case, the micro-fracture cloud is characterized by a Hurst exponent of 0.4 for DEM joints S5 to S8 and $H \approx 0.8$, for DEM joints S1 to S4, closely resembling the roughness exponent measured by Ponson et al. [⁴²] and Morel et al. [⁴³].

Finally, Figure 11 shows the evolution of micro-fracture density as function of the perpendicular distance ($\pm z$ direction) to the average joint surface. Note that the fracture density logarithmically decreases with the distance from the average joint plane. This agrees with the experimental results reported by Vermilye et al. [⁴⁴].
<figure><img src="image_15.png" /><figcaption>**Figure 10** Evolution of micro-fracture density perpendicular to the mean plane of the fracture (±z direction). Each color represents a DEM joint.</figcaption></figure>

**5. Conclusions**

In summary, our study underscores the existence of spatial-temporal correlations among micro-fractures in self affine rough joints subjected to shear.

* Our approach is based upon a DEM model that faithfully reproduces mechanical behaviors akin to real-world scenarios, encompassing compression, tension, and shear responses. The DEM-generated joints underwent shear under constant normal load conditions at a consistent shear speed.

* To explore the temporal relationships between fractures, we computed the waiting times (Δt) separating successive fractures. We distinguished between two critical cases: 1) interacting micro-cracks, occurring in proximity within a confined region, and 2) non-interacting micro-cracks, where fractures emerged across the entire joint volume. The probability of encountering two fractures within a specified time lapse follows a power-law distribution akin to Omori's law, commonly employed in earthquake studies. The exponent for interacting micro-fractures ranged between 0.716 and 0.869, while for all micro-fractures, it spanned from 1.458 to 1.895. Notably, we observed a dependence of these temporal correlations on the roughness characteristics of the joints for all micro-fracturing events, this dependency is more visible between the smoothest (DEMS4) and roughest (DEMS5) joints. The events corresponding to interacting cracks depend only slightly on the joint roughness. They are too close to each other to feel the effect of roughness.

* Observations reveals also a propensity for micro-fractures to cluster during shearing. The probability of finding two fractures separated by a given distance (Δd) follows a power law of this distance with an exponent q. This exponent fluctuates from 0.459 for the roughest joint (DEM S5) to 0.843 for the smoothest joint (DEM S4). Micro-fracture clouds exhibit fractal characteristics and the hurst exponent of the fractal set of micro-fracture is equivalent to exponent q.
* It is worth highlighting that our findings also show remarkable agreement with the experimentally reported values in the literature for the *p* and *q* exponents.

**6. Acknowledgements**

The authors are grateful with CONACYT for the financial support. Financial support from the investissement d'avenir French programme (ANR-16-IDEX-0002) to G. Pijaudier-Cabot through the E2S hub Newpores is gratefully acknowledged.

**7. References**

1. Byerlee J. Friction of rocks. *Pure Appl Geophys PAGEOPH*. 1978;116(4-5):615-626.
doi:10.1007/BF00876528

2. Swenson D, Hardeman B. The effects of thermal deformation on flow in a jointed geothermal reservoir. *Int J rock Mech Min Sci Geomech Abstr*. 1997;34(3-4):445. doi:10.1016/S1365-1609(97)00285-2

3. Pandey SN, Chaudhuri A, Kelkar S. A coupled thermo-hydro-mechanical modeling of fracture aperture alteration and reservoir deformation during heat extraction from a geothermal reservoir. *Geothermics*. 2017;65:17-31. doi:10.1016/j.geothermics.2016.08.006

4. Barton N, Bandis S, Bakhtar K. Strength, deformation and conductivity coupling of rock joints. *Int J Rock Mech Min Sci*. 1985;22(3):121-140. doi:10.1016/0148-9062(85)93227-9

5. Barton N, de Quadros EF. Joint aperture and roughness in the prediction of flow and groutability of rock masses. *Int J Rock Mech Min Sci*. 1997;34(3-4):252.e1-252.e14. doi:10.1016/S1365-1609(97)00081-6

6. Barton N, Choubey V. The shear strength of rock joints in theory and practice. *Rock Mech Felsmechanik Mecanique des Roches*. 1977;10(1-2):1-54. doi:10.1007/BF01261801

7. Bandis S, Lumsden AC, Barton NR. Experimental Studies of Scale Effects on the Shear Behaviour of Rock Joints. *Rock Mech Min Sci*. 1981;18:1-21. doi:10.1016/0148-9062(81)90262-x

8. Grasselli G, Egger P. Constitutive law for the shear strength of rock joints based on three-dimensional surface parameters. *Int J Rock Mech Min Sci*. 2003;40(1):25-40. doi:10.1016/S1365-1609(02)00101-6

9. Bandis SC, Lumsden AC, Barton NR. Fundamentals of rock joint deformation. *Int J Rock Mech Min Sci*. 1983;20(6):249-268. doi:10.1016/0148-9062(83)90595-8

10. Renard F, Candela T, Bouchaud E. Constant dimensionality of fault roughness from the scale of micro-fractures to the scale of continents. *Geophys Res Lett*. 2013;40(1):83-87.
doi:10.1029/2012GL054143

11. Kulatilake PHSW, Balasingam P, Park J, Morgan R. Natural rock joint roughness quantification through fractal techniques. *Geotech Geol Eng*. 2006;24(5):1181-1202. doi:10.1007/s10706-005-1219-6

12. Ponson L, Auradou H, Pessel M, Lazarus V, Hulin JP. Failure mechanisms and surface roughness
statistics of fractured Fontainebleau sandstone. *Phys Rev E - Stat Nonlinear, Soft Matter Phys.* 2007;76(3):1-7. doi:10.1103/PhysRevE.76.036108

13. Mourot G, Morel S, Bouchaud E, Valentin G. Scaling properties of mortar fracture surfaces. *Int J Fract.* 2006;140(1-4):39-54. doi:10.1007/s10704-005-3471-4

14. Candela T, Renard F, Bouchon M, et al. Characterization of Fault Roughness at Various Scales: Implications of Three-Dimensional High Resolution Topography Measurements. *Pure Appl Geophys.* 2009;166(10-11):1817-1851. doi:10.1007/s00024-009-0521-2

15. Neuville a., Toussaint R, Schmittbuhl J. Fracture aperture reconstruction and determination of hydrological properties: A case study at Draix (French Alps). *Hydrol Process.* 2012;26(14):2095-2105. doi:10.1002/hyp.7985

16. Neuville A, Toussaint R, Schmittbuhl J. Hydrothermal coupling in a self-affine rough fracture. *Phys Rev E.* 2010;82(3):036317. doi:10.1103/PhysRevE.82.036317

17. Schmittbuhl J, Vilotte J-P, Roux S. Dynamic friction of self-affine surfaces. *J Phys II Fr.* 1994;4:225-237.

18. Neuville A, Toussaint R, Schmittbuhl J. Hydraulic transmissivity and heat exchange efficiency of open fractures: A model based on lowpass filtered apertures. *Geophys J Int.* 2011;186(3):1064-1072. doi:10.1111/j.1365-246X.2011.05126.x

19. Varela Valdez A, Morel S, Marache A, Hinojosa M, Riss J. Influence of fracture roughness and micro-fracturing on the mechanical response of rock joints: a discrete element approach. *Int J Fract.* 2018;213(2):87-105. doi:10.1007/s10704-018-0308-5

20. Yang Z-Y, Taghichian A, Huang G-D. On the applicability of self-affinity concept in scale of three-dimensional rock joints. *Int J Rock Mech Min Sci.* 2011;48(7):1173-1187. doi:10.1016/j.ijrmms.2011.06.010

21. Wei M-Y, Liu H-H, Li L-C, Wang E-Y. A fractal-based model for fracture deformation under shearing and compression. *Rock Mech Rock Eng.* 2013:1539-1549. doi:10.1007/s00603-013-0367-x

22. Lo TS, Koplik J. Channeling and stress during fluid and suspension flow in self-affine fractures. *Phys Rev E - Stat Nonlinear, Soft Matter Phys.* 2014;89(2):1-14. doi:10.1103/PhysRevE.89.023010

23. Giacomini A, Buzzi O, Ferrero AM, Migliazza M, Giani GP. Numerical study of flow anisotropy within a single natural rock joint. *Int J Rock Mech Min Sci.* 2008;45(1):47-58. doi:10.1016/j.ijrmms.2007.04.007

24. Yabe Y. Evolution of source characteristics of AE events during frictional sliding. *Earth, Planets Sp.* 2008;60(12):8-11.

25. Goebel THW, Candela T, Sammis CG, Becker TW, Dresen G, Schorlemmer D. Seismic event distributions and off-fault damage during frictional sliding of saw-cut surfaces with pre-defined roughness. *Geophys J Int.* 2014;196(1):612-625. doi:10.1093/gji/ggt401

26. Meng F, Zhou H, Wang Z, et al. Influences of Shear History and Infilling on the Mechanical Characteristics and Acoustic Emissions of Joints. *Rock Mech Rock Eng.* 2017;50(8):2039-2057. doi:10.1007/s00603-017-1207-1
27. Moradian ZA, Ballivy G, Rivard P, Gravel C, Rousseau B. Evaluating damage during shear tests of rock joints using acoustic emissions. *Int J Rock Mech Min Sci.* 2010;47(4):590-598. doi:10.1016/j.ijrmms.2010.01.004

28. Rudajev V, Vilhelm J, Kozak J, Lokajicek T. Statistical precursors of instability of loaded rock samples based on acoustic emission. *Int J Rock Mech Min Sci Geomech Abstr.* 1996;33(7):478-743. doi:http://dx.doi.org/10.1016/0148-9062(96)00023-X

29. Utsu T, Ogata Y, Matsu'ura RS. The centenary of the Omori Formula for a decay law of aftershock activity. *J Phys Earth*. 1995;43:1-33.

30. Bunger AP, Kear J, Dyskin A V., Pasternak E. Sustained acoustic emissions following tensile crack propagation in a crystalline rock. *Int J Fract.* 2015;193(1):87-98. doi:10.1007/s10704-015-0020-7

31. Ojala IO, Main IG, Ngwenya BT. Strain rate and temperature dependence of Omori law scaling constants of AE data: Implications for earthquake foreshock-aftershock sequences. *Geophys Res Lett.* 2004;31(24):1-5. doi:10.1029/2004GL020781

32. Smirnov VB, Ponomarev A V., Stanchits SA, et al. Laboratory Modeling of Aftershock Sequences: Stress Dependences of the Omori and Gutenberg-Richter Parameters. *Izv Phys Solid Earth*. 2019;55(1):124-137. doi:10.1134/S1069351319010105

33. Kagan YY, Knopoff L. Stochastic synthesis of earthquake catalogs. *J Geophys Res Solid Earth*. 1981;86(B4):2853-2862. doi:https://doi.org/10.1029/JB086iB04p02853

34. Itasca Inc. PFC3D. 2008.

35. Niccolini G, Carpinteri A, Lacidogna G, Manvello A. Acoustic emission monitoring of the Syracuse Athena temple: Scale invariance in the timing of ruptures. *Phys Rev Lett.* 2011;106(10). doi:10.1103/PhysRevLett.106.108503

36. Schmid A, Grasso JR. Omori law for eruption foreshocks and aftershocks. *J Geophys Res Solid Earth*. 2012;117(7). doi:10.1029/2011JB008975

37. Aguirre D, Hinojosa JAM. Microestructura, propagación y emisión acústica de grietas en dos papeles de celulosa. *Rev Mex Fis.* 2015;61(August):254-260.

38. Karabulut H, Schmittbuhl J, Özalaybey S, et al. Evolution of the seismicity in the eastern Marmara Sea a decade before and after the 17 August 1999 Izmit earthquake. *Tectonophysics*. 2011;510(1-2):17-27. doi:10.1016/j.tecto.2011.07.009

39. Grassberger P, Procaccia I. Measuring the strangeness of strange attractors. *Phys D Nonlinear Phenom*. 1983;9(1-2):189-208. doi:10.1016/0167-2789(83)90298-1

40. Godano C, Tosi P, Derubeis V, Augliera P. Scaling properties of the spatio-temporal distribution of earthquakes: A multifractal approach applied to a Californian catalogue. *Geophys J Int.* 1999;136(1):99-108. doi:10.1046/j.1365-246X.1999.00697.x

41. Lefort V, Pijaudier-Cabot G, Grégoire D. Analysis by Ripley's function of the correlations involved during failure in quasi-brittle materials: Experimental and numerical investigations at the mesoscale. *Eng Fract Mech.* 2015;147:449-467. doi:10.1016/j.engfracmech.2015.07.046

42. Ponson L, Auradou H, Pessel M, Lazarus V, Hulin JP. Failure mechanisms and surface roughness statistics of fractured Fontainebleau sandstone. *Phys Rev E*. 2007;(76):1-7.
doi:10.1103/PhysRevE.76.036108

43. Morel S, Bonamy D, Ponson L, Bouchaud E. Transient damage spreading and anomalous scaling in mortar crack surfaces. *Phys Rev E*. 2008;78(1):016112. doi:10.1103/PhysRevE.78.016112

44. Vermilye JM, Scholz CH. The process zone: A microstructural view of fault growth. *J Geophys Res Solid Earth*. 1998;103(B6):12223-12237. doi:10.1029/98JB00957