# Role of heterogeneities in jointing (fracturing) of geological media: Numerical analysis of fracture mechanisms

Alexandre Chemenda, Julien Ambre, Junchao Chen, Jinyang Fan, Jean-Pierre Petit, Deyi Jiang

► To cite this version:

Alexandre Chemenda, Julien Ambre, Junchao Chen, Jinyang Fan, Jean-Pierre Petit, et al.. Role of heterogeneities in jointing (fracturing) of geological media: Numerical analysis of fracture mechanisms. Comptes Rendus Géoscience, 2018, 350 (8), pp.452-463. 10.1016/j.crte.2018.06.009. hal-01983594

HAL Id: hal-01983594

https://hal.science/hal-01983594

Submitted on 28 Aug 2020

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
Tectonics, Tectonophysics

# Role of heterogeneities in jointing (fracturing) of geological media: Numerical analysis of fracture mechanisms

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

Check for
updates

Alexandre I. Chemenda<sup>a,*</sup>, Julien Ambre<sup>a</sup>, Junchao Chen<sup>a,b</sup>, Jinyang Fan<sup>b</sup>, Jean-Pierre Petit<sup>c</sup>, Deyi Jiang<sup>b</sup>

## ARTICLE INFO

**Article history:**

Received 29 June 2018

Accepted after revision 29 June 2018

Available online 6 September 2018

Handled by Vincent Courtillot

**Keywords:**

Fracture

Damage

Joints

Finite-difference modeling

Plumose

Plasticity theory

AMS codes: 65, 74, 86

## ABSTRACT

Fracture process is investigated using finite-difference simulations with a new constitutive model. It is shown that both geometry and fracture mechanism itself depend on the preexisting heterogeneities that are stress concentrators. In the brittle regime (low pressure, P), Mode-I fractures propagate normal to the least stress $\sigma_3$ from the imposed weak zones. At high P, shear deformation bands are formed oblique to $\sigma_3$. At intermediate values of P, the fracture process involves both shear banding and tensile cracking and results in the initiation and propagation of pure dilation bands. The propagating band tip undulates, reacting on the failure mechanism changes, but its global orientation is normal to $\sigma_3$. The $\sigma_3$-normal fractures are joints. There are thus two types of joints resulting from Mode-I cracking and dilation banding, respectively. The obtained numerical results are in good agreement with and explain the results from previous similar experimental study. © 2018 Published by Elsevier Masson SAS on behalf of Académie des sciences. This is an open access article under the CC BY-NC-ND license ([http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/)).

## 1. Introduction

There are two principal types of fractures in geological media and notably in sedimentary basins (reservoirs), joints and shear fractures. The principal difference between them is that shear fractures accommodate shear displacements, whereas along joints only normal to joints displacements are possible. Another difference is in the orientation of fractures compared to the acting (during fracturing) stress directions. Joints are parallel to the most compressive stress $\sigma_1$ and normal to the least compressive or the most tensile stress $\sigma_3$. Joints thus make angle $\psi = 0^\circ$ with $\sigma_1$ and typically are organized in regular parallel

(tabular) sets or orthogonal networks. Shear fractures are oblique to both $\sigma_1$ and $\sigma_3$ and typically form conjugate networks, oriented to $\sigma_1$ at angles $\psi = \pm 10^\circ - 30^\circ$. All fracture types are parallel to the intermediate principal stress $\sigma_2$.

### 1.1. Theoretical background of fracture formation

The mechanism of fracture formation remains not completely clear. This is particularly true for joints, which is reflected in different names used for this fracture type: tensile, extension, opening mode or cleavage fractures. Most commonly, they are believed to form (or more exactly, to propagate from a preexisting flaw) according to Mode-I mechanism (Fig. 1) stemming from Linear Elastic Fracture Mechanics (LEFM), (e.g., Lawn, 1993). Accordingly, shear fractures are believed to be Mode-II fractures (Fig. 1). From a theoretical point of view, LEFM is applicable

<footnote><sup>a</sup> University of Côte d'Azur, CNRS, OCA, IRD, Géoazur, 250, rue Albert-Einstein, Sophia Antipolis, 06560 Valbonne, France</footnote>

<footnote><sup>b</sup> State Key Laboratory for the Coal Mine Disaster Dynamics and Controls, Chongqing University, 400044 Chongqing, China</footnote>

<footnote><sup>c</sup> University of Montpellier, Géosciences Montpellier, place Eugène-Bataillon, 34095 Montpellier cedex 05, France</footnote>

<footnote>* Corresponding author.
E-mail address: chem@geoazur.unice.fr (A.I. Chemenda).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Fracture modes.</figcaption></figure>

only when the size of the near tip zone (process zone) undergoing inelastic (irreversible) deformation or damage is small compared to other geometric dimensions of the problem (fracture length, size of the deforming object). This assumption can hold only at a relatively low pressure P (or stress level), when deformation and fracture occur in a purely brittle regime. However, this is certainly not the case in the experimental tests at a relatively high P, when damage (revealed by acoustic emissions) involves practically the entire deforming specimen prior to the fracture (e.g., Fortin et al., 2006). In this “ductile” regime, the failure occurs via localization of the initially distributed damage (inelastic deformation) within a narrow band that then can become a fracture. The theoretical framework for addressing the mechanical instability resulting in deformation localization constitutes the bifurcation theory applied to elasto-plasticity models (Chemenda, 2007; Issen and Rudnicki, 2000; Rudnicki and Rice, 1975; Vardoulakis and Sulem, 1995). This theory is completely different from LEFM (which does not explicitly take into account the inelastic deformation), and explains very well the formation of various types of shear deformation bands at the origin of shear fractures. However, the prediction of dilation banding (which can be at the origin of joints, see below) does no go straightforward from this theory when realistic material parameters (notably the dilatancy factor) are adopted.

## 1.2. Insights from experimental studies

The transition from very brittle to “ductile” fracturing with increasing pressure has been observed in numerous experimental tests (e.g., Bésuelle et al., 2000; Fortin et al., 2006; Wong et al., 1997) including the conventional (axisymmetric) extension tests on specimens made of Granular Rock Analogue Material (GRAM), which behaves as real rocks, but is much weaker than hard rocks (Nguyen et al., 2011). The dog-bone and cylindrical specimens are first loaded hydrostatically in these tests to a given confining pressure P that is equal to σ₁ and σ₂, and then unloaded in the axial direction

that becomes σ₃, P being maintained constant. With increasing P, specimen fracture occurs at progressively smaller |σ₃| (which is negative), but the failure plane orientation does not change until σ₃ reaches approximately zero (Fig. 2a, b): this plane is parallel to σ₁ (ψ ≈ 0°), Fig. 2b, which corresponds to jointing. With further P increase, ψ increases (Fig. 2b), which corresponds to shear fracturing.

Although the orientation of failure planes in the stress space in GRAM experiments is the same (ψ ≈ 0°) within the P range of joints formation, the “ductility” of the failure process clearly increases with P, which is expressed in the morphology of the forming fracture surfaces (Chemenda et al., 2011a). At low P, these surfaces are rather smooth, while at high P, they are rough, with the rugosities being organized in nice plumose features (Fig. 2c) so typical for natural joints (e.g., Bahat, 1991; Pollard and Aydin, 1988; Savalli and Engelder, 2005; Woodworth, 1896). The damage/failure process occurs under these pressures within a band of a finite thickness comparable to the plumose relief. This process is therefore rather ductile and is not limited to the opening of a “zero”-thickness plane, which would be the case during Mode-I fracturing. SEM observations of the failure zones forming at high P in GRAM specimens reveal a several grains-thick irregular wavy band of damaged material with increased porosity (Chemenda et al., 2011a). This σ₃-normal band can be called a pure dilation band. The conclusion made by Chemenda et al. (2011a, b) is that the jointing mechanism changes from Mode-I cracking to dilatancy/dilation banding with increasing P. With further pressure increase, failure occurs via the formation of shear-dilatant deformation bands.

The dilation banding clearly involves ductile, but also tensile damage/failure mechanisms. This is the reason why the dilatancy banding could not be correctly predicted theoretically from deformation bifurcation analysis within an elasto-plasticity framework as mentioned above. In this paper, we investigate this process numerically by considering both brittle-tensile (BT) and shear-ductile (SD) failure/damage mechanisms.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Results from conventional (axisymmetric) extension tests of Granular Rock Analogue Material, GRAM (from Chemenda et al., 2011a). (a) Minimum nominal stress σ₃ at fracturing/rupture vs confining pressure or σ₁. (b) Angle ψ between the rupture plane and σ₁ direction vs σ₁ at failure. (c) Surface of the fracture (opened rupture zone) formed at σ₁ = σ₂ = 0.6 MPa. The initiation point (zone) of the plumose morphology is due to the stress concentration caused by the rigid support of the LVDT gage attached to the specimen (to the jacket) at that location.</figcaption></figure>

## 2. Constitutive formulation

The adequate constitutive formulation is the most important ingredient in theoretical analysis and numerical modeling of inelastic deformation (damage) and failure of geomaterials. Rock mechanics experimental studies with different loading configurations (considerably accelerated recently (e.g., Ingraham et al., 2013; Li et al., 2018; Ma et al., 2017), clearly show the importance of the loading type (controlled by the Lode angle), which affects the mechanical response of rocks. A new constitutive model has been proposed based on these experimental data (Chemenda and Mas, 2016). The model includes a three-invariant yield function

$$F = \bar{\tau} - \bar{\tau}_e(\sigma) \left\{ \frac{1}{2} \left[ 1 + \cos(3\theta) + R(\sigma)^{\frac{1}{n}}(1 - \cos(3\theta)) \right] \right\}^n \quad (1)$$

and the related plastic potential (n is the constant scalar that can take values from -0.1 to -1). The bifurcation analysis within the framework of this constitutive model led to the predictions of rock failure regimes that are in good agreement with the experimental data, and that are different from the predictions made based on the classical 2-invariant theories like that of Drucker-Prager (Chemenda and Mas, 2016). Here $\bar{\tau}$ is von Mises's stress; $R = \bar{\tau}_c/\bar{\tau}_e$; $\bar{\tau}_c$ and $\bar{\tau}_e$ are functions of the mean stress σ defined from the conventional compression and extension tests, respectively; θ is the Lode angle. This constitutive model is rather complicated and has not been implemented into a numerical code as yet. For this reason, we use here much simpler, but also 3-invariant Coulomb-type theory. The yield function is taken in the form

$$F = A\sigma_1 - \mu_c\sigma_3 - k_c \quad (2)$$

It does not include σ₂, which is a simplification of reality that can be justified by the experimental fact that the rock properties depend much stronger on σ₁ and σ₃ than on σ₂ (e.g., Ingraham et al., 2013; Ma et al., 2017). The plastic potential function reads

$$G = A\sigma_1 - \mu_c\sigma_3 \quad (3)$$

where $\mu_c$, $k_c$, $\mu_β$ are the friction, cohesion and dilatancy parameters (not coefficients); A takes value 0 or 1 for the brittle-tensile (BT) and shear ductile (SD) failure/damage mechanisms, respectively. The initial values of $\mu_c$ and $k_c$ as well as of the tensile strength $\sigma_t$ ($\mu_c^0, k_c^0$ and $\sigma_t^0$) are defined from the GRAM experimental data (from the simplified $\sigma_3(\sigma_1)$ curve in Fig. 3a) to be

$$\begin{aligned} \mu_c^0 &= 1, A = 0, k_c = -\sigma_t & \text{if } \sigma_3 &= -\sigma_t \\ \mu_c^0 &= 2, A = 1, k_c = 0.65 \text{ MPa} & \text{if } \sigma_3 &>& -\sigma_t, \end{aligned} \quad (4)$$

The much higher brittleness of the BT (than that of the SD) mechanism is taken into account in the following softening rules:

$$\begin{aligned} \sigma_t &= \sigma_t^0 \left(1 - \frac{\varepsilon_t}{\varepsilon_0}\right); & k_c &= k_c^0 \left(1 - \frac{\varepsilon}{\varepsilon_0}\right); \\ \mu_c &= \mu_c^0 + \left(\mu_c^{\text{end}} - \mu_c^0\right) \frac{\varepsilon}{\varepsilon_0} \end{aligned} \quad (5)$$

where $\sigma_t^0 = 0.07 \text{ MPa}$; $\varepsilon_t^0 = 0.5 \times 10^{-3}$ and $\varepsilon^0 = 1 \times 10^{-2}$ ($\varepsilon_t^0 \ll \varepsilon^0$); $\mu_c^{\text{end}} = 1.7$, $\varepsilon_t = \int d\varepsilon_3^p$ and $\varepsilon_t = \varepsilon_t + \int d\varepsilon_1^p$ (the superscript “p” stands for plastic or inelastic). Note that $\sigma_1$ in (5) evolves only with the accumulated inelastic extension strain, which is a sum of strains resulting from BT and SD deformations. This is the way how both failure mechanisms are coupled. According to this formulation, the tensile strength decreases even during SD deformation along with the cohesion decrease, meaning that these two parameters are related, both characterizing the pressure-independent part of the strength. According to the results of the processing of a large experimental data set (Mas and Chemenda, 2015), the dilatancy parameter $\beta_c$ evolves (reduces) with $\varepsilon$ in the same way as $\mu_c$, see Eq. (5).

The elastic properties of the model both before and after reaching the yield condition are described by Hooke's equations with Poisson's ratio $\nu = 0.25$ and Young's modulus $E = 6.7 \times 10^8 \text{ Pa}$, characterizing GRAM material (Nguyen et al., 2011).

The formulated constitutive model has been implemented into the finite-difference dynamic time-matching explicit code Flac3D (Itasca, 2012) used for numerical modeling in this work. Numerous simulations were carried out in quasi-static, large-strain mode with different model parameters. Below we report representative results.
<figure><img src="image_5.png" /><figcaption>Fig. 3. (a) Initial yield surface/envelope or function used in the numerical simulations (simplified plot in Fig. 2a). Points 1 to 4 correspond to the stress states applied in the numerical models. (b) Setup of numerical modeling. The models were pre-stressed to the initial stresses $σ_1$ and $σ_2$ ($σ_2 = σ_1$) corresponding to points 1 to 4 in Fig. 3a, and to a sufficiently large $σ_3$ for the model to remain in a purely elastic state. Then $σ_3$ was progressively decreased (by applying velocities $V_3$ to the $σ_3$-normal boundaries), while maintaining the initial values of $σ_1$ and $σ_2$ by applying velocities $V_1$ and $V_2$ as shown in this figure. The model size is 5 × 10 cm (200 × 400 numerical zones) in the $σ_1$ and $σ_3$ directions, and one numerical zone, in $σ_2$ direction. To take into account in the models heterogeneities typical for real rocks, two scale property heterogeneities were introduced in different models (indicated in the captions of the figures showing the corresponding numerical results): (1) Gauss distribution in the numerical zones of the Poisson’s ratio with a mean value of 0.25 and standard deviation of 0.02; (2) a weak zone of a larger scale (including several numerical zones) shown in this figure. The strength in this zone is indicated in the captions of the figures presenting the numerical results.</figcaption></figure>

### 3. Setup of numerical modeling

The modeling setup (Fig. 3b) is two-dimensional in the sense that the model thickness in the $σ_2$-direction is only one numerical zone (or 0.25 mm), but the strain rate in this direction is not zero (see below). The model is first pre-stressed to the initial stresses $σ_1 = σ_2$ and to a sufficiently large $σ_3$ for the stress state to be above the initial yield envelope in Fig. 3a (this corresponds to the purely elastic mechanical state). Then, $σ_3$ is progressively decreased by applying velocities $V_3$ to the $σ_3$-normal model boundaries, while maintaining the initial values of $σ_1$ and $σ_2$ by applying velocities $V_1$ and $V_2$ as shown in Fig. 3b. In reality, maintained and monitored are averages over each boundary or nominal stresses as it is typically done in the experimental tests. Since $V_2$ is not zero, the models are not plane-strain.

### 4. Results

We report here the results corresponding to the four points shown on the yield envelope $σ_3(σ_1)$ in Fig. 3a. The first point (point 1, $σ_1$ = 0.3 MPa; $σ_3$ = $-σ_1^0$ = $-0.07$ MPa) is located in the middle of the BT segment of the envelope. Point 2 ($σ_1$ = 0.54 MPa; $σ_3$ = $-0.055$ MPa) is located on the SD segment close to the transition to the BT segment. Point 3 ($σ_1$ = 0.6 MPa; $σ_3$ = $-0.025$ MPa) is on the SD segment close to $σ_3$ = 0 ($σ_3$ < 0), and point 4 ($σ_1$ = 0.66 MPa; $σ_3$ = 5 × 10³ Pa) is also on the SD segment with $σ_3$ > 0.

#### 4.1. Point 1 in Fig. 3a, $σ_1 = σ_2$ = 0.3 MPa

##### 4.1.1. Homogeneous model

In this case, the tensile failure occurs exactly at (along) one of the $σ_3$-normal ends of the model due to the end effect. Although very small, this effect exists in the

numerical models, and occurs where the extension is applied to the model, i.e. at the $σ_3$-normal boundaries.

##### 4.1.2. "Homogeneously heterogeneous" model with small variations of the Poisson's ratio $v$

The Gauss distribution of $v$ is applied in this model with a mean value of 0.25 and a standard deviation of 0.02. This leads to a non-uniform stress distribution during the loading, which eliminates the boundary (end) affects because the stress perturbations within the model are larger than those (very small) caused by these effects.

Fig. 4 shows the evolution of fracturing, which occurs exclusively according to the BT mechanism as expected (Fig. 4e); the SD mechanism was not activated at all. The formed tensile fractures constitute a generally normal to $σ_3$ but rather irregular set. This is because the model is almost homogeneous and therefore fracturing is initiated simultaneously in different points resulting in stress drops and perturbations of the stress field. Therefore, the initiated fractures propagate in the heterogeneous field, which explains their irregular geometry. To avoid such multiple stress perturbations, in the following model a small weak zone is introduced.

##### 4.1.3. Model with weak zone

This zone is shown in Fig. 3b and can be recognized in Fig. 5a. The strength parameters in this zone are reduced, as indicated in the caption of Fig. 5; the other parameters are the same as in the rest of the model. One can see that the fracturing occurs in this case in a completely different way compared to the previous model in Fig. 4. There is only one fracture that is initiated at the weak zone and propagates strictly normal to $σ_3$. Propagation occurs due to the tensile stress concentration at the fracture tip, where it reaches $σ_t$ (Fig. 6b). This Mode-I-type propagation occurs exclusively due to the BT failure mechanism,
<figure><img src="image_6.png" /><figcaption>Fig. 4. (a–d) Evolution of the inelastic strain, $\varepsilon^p$, pattern during vertical unloading ($\sigma_3$ reduction) of the model at constant $\sigma_2 = \sigma_1 = 0.3$ MPa, corresponding to point 1 in Fig. 3a. Poisson’s ratio $\nu$ varies in the numerical zones according to Gauss distribution. (e) Mechanical state (elastic and elastic-plastic), damage/ failure or inelastic deformation mechanism (tensile or shear) and its history for the deformation stage in (d). Different colors in this image correspond to different damage/failure mechanisms and their succession. For example, “tensile-n (n means now), shear-p, tensile-p (p means past)” indicates that the corresponding zones undergo tensile damage at the stage shown, but they already underwent shear and tensile damage during the previous deformation stages (in the past). Zones indicated as “elastic” did not undergo inelastic straining at all.</figcaption></figure>

SD mechanism being operated only within the weak zone (Figs. 6c, d).

## 4.2. Point 2 in Fig. 3a, $\sigma_1 = \sigma_2 = 0.54$ MPa

### 4.2.1. Homogeneous model

A dense network of conjugate shear bands is first initiated throughout the entire model (Fig. 7a). Then, most initial bands die, while the remaining ones accommodate

growing inelastic deformation. This process results in a few principal (well-developed) shear bands (Fig. 7c). One of them then is transferred into three-segment fracture when the BD mechanism comes in play (Fig. 7d, e).

### 4.2.2. Homogeneously heterogeneous model

Before fracturing, this model undergoes distributed BT, but mostly SD distributed damage (inelastic deformation), Fig. 8a. Then, an oblique fracture is initiated and

<figure><img src="image_7.png" /><figcaption>Fig. 5. Evolution of the $\varepsilon^p$ pattern in the model with Gaussian $\nu$ distribution and a weak zone with the reduced initial strength parameters: $\sigma_t = 0.03$ MPa; $k_c = 0.1$ MPa. The initial stresses and boundary conditions are the same as in the previous model.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. The inelastic strain $\dot{\epsilon}^p$ (a), the least stress $\sigma_3$ (b), and the mechanical state (c) patterns at the same evolutionary stage of the model in Fig. 5; (d) corresponds to the last stage of fracture propagation in Fig. 5e.</figcaption></figure>

propagates due to both BT and SD mechanisms (Fig. 8b-d). The oblique fracturing process is followed by an almost $\sigma_3$-normal one (Fig. 7d).

### 4.2.3. Model with weak zone

The introduction of the weak zone again changes radically the fracture process and geometry (Fig. 9). The fracture propagating from the weak zone does not follow exactly a linear trajectory in this case. Its tip slightly undulates during propagation, resulting in globally $\sigma_3$-normal macro fracture orientation, whose formation involves both BT and SD mechanisms (Fig. 9). One can recognize both the process (blue dots) and damage (green dots) zones on the maps of damage distribution in Fig. 9 (the second row in this figure).

### 4.3. Point 3 in Fig. 3a, $\sigma_1 = \sigma_2 = 0.6$ MPa

#### 4.3.1. Homogeneous model

The deformation pattern in this case is similar to that in Fig. 7, which is explained by the fact that, in both cases, the imposed stress states correspond to the SD domain.

#### 4.3.2. Model with weak zone (Fig. 10)

The result is similar to that in Fig. 9, but the process and damage zones related to the fracture (or more exactly to the deformation band) propagation are much larger. They develop exclusively through an SD mechanism and involve almost the entire model (the last row in Fig. 10). At the band/fracture tip, $\sigma_3$ is tensile. Its maximum magnitude reaches around 0.06 MPa, which is somewhat less than $\sigma_t^0 = 0.07$ MPa. The nominal stress (which would be

measured at the model (sample) end in the experimental test) is $\sigma_3 = -0.0011$ MPa, which is much less in magnitude than $\sigma_t^0$.

Within the band, both BT and SD mechanisms are active. As in the previous case (Fig. 9), the fracture/band tip undulates during propagation. This time, the reason for this is clear from the octahedral strain rate $\dot{\epsilon}$ and $\sigma_t$ patterns in Fig. 10. During generally $\sigma_3$-normal band tip propagation, the initiation of oblique shear bands occurs. The fracture tip tends to follow these (oblique) directions alternately, keeping general trajectory normal to $\sigma_3$. In the next model, it is the oblique direction that is definitively followed by the deformation band.

### 4.4. Point 4 in Fig. 3b, $\sigma_1 = \sigma_2 = 0.66$ MPa

#### 4.4.1. Model with weak zone

The same weak zone as in the previous model leads in this case to the initiation of the $\sigma_3$-normal extension band (that can hardly be recognized in Fig. 11) and of two conjugate shear bands (Fig. 11a), the latter propagating much faster than the former. Then one shear band dies and the other one cuts the entire model and becomes a shear fracture (Fig. 11e). One can see a large process zone in front of the propagating band (the second row in Figs. 11b-d).

## 5. Concluding discussion

The reported results show that not only the geometrical attributes of fracture in geomaterials, but the fracture mechanism itself strongly depend on the presence of heterogeneities. Joints, which are believed to be Mode-I
<figure><img src="image_9.png" /><figcaption>Fig. 7. Evolution of the $\varepsilon^p$ and mechanical state patterns in a homogeneous model corresponding to point 2 ($\sigma_2 = \sigma_1 = 0.54$ MPa) in Fig. 3a.</figcaption></figure>

fractures, are supposed to form exclusively when the tensile stress $\sigma_3$ reaches the tensile strength $\sigma_t$ (or more exactly, when the stress intensity factor reaches a critical value). This concept works fairly well for purely tensile fractures forming according to the brittle-tensile (BT) deformation/fracture mechanism. This mechanism is possible at a relatively low far-field stress level (or pressure $P$), when the inelastic process zone near the fracture tip is very small as in Figs. 5 and 6. However, when $P$ is high enough (compared to $\sigma_t^0$), joints can form at $|\sigma_3| < <\sigma_t^0$ (Figs. 8 and 9). Here $\sigma_3$ is not the stress at the micro- (granular) scale, but the “far-field” stress applied to the boundaries of the object (specimen or structure). This corresponds to the nominal stress $\sigma_3$, which can be

measured experimentally in the laboratory. At the micro-scale, $|\sigma_3|$ is tensile, but lower than $\sigma_t^0$. Purely tensile, $\sigma_3$-normal fractures cannot propagate in this case. The damage (inelastic deformation) involves a wide zone (almost the entire model) emanating from the fracture front (Figs. 10b–d). The failure of the material occurs in this zone according to a shear ductile (SD) mechanism, leading to the initiation of oblique shear bands. The fracture tip tends to follow these oblique directions alternately, but keeps general $\sigma_3$-normal orientation because of the involvement of the BT mechanism in the failure process. As a result, the fracture zone is not linear. It is thicker and more heterogeneous than at lower $P$. This zone corresponds to the dilation bands (embryonic joints) obtained
<figure><img src="image_10.png" /></figure>

Mechanical state, damage/failure mechanism and its history

<figure><img src="image_11.png" /><figcaption>Fig. 8. Evolution of the $\varepsilon^p$ and mechanical state patterns in the model with a Gaussian distribution of $v$. The initial stresses and boundary conditions are the same as in the previous model (the stress state corresponds to point 2 in Fig. 3a).</figcaption></figure>

experimentally in Chemenda et al. (2011a). The walls of the opened band cannot be smooth (as in the case of the “ideal” Mode-I fracture) and should bear traces (irregularities) reflecting complex brittle-ductile failure process.

These traces are likely at the origin of a plumose morphology of the walls of natural (as well as experimental, Fig. 2c) joints, which is usually interpreted as a result of mixed mode (I and III, Fig. 1) fracturing in the sense of
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Evolution of the $\varepsilon^p$ and mechanical state patterns in the model with a Gaussian distribution of $\nu$ and a weak zone with following strength parameters: $\sigma_t = 0.03$ MPa; $k_c = 0.3$ MPa. The initial stresses and boundary conditions are the same as in the two previous models (the stress state corresponds to point 2 in Fig. 3a).</figcaption></figure>

Linear Elastic Fracture Mechanics, LEFM, (Pollard and Aydin, 1988). Both opening and shear are indeed present in the described dilatancy banding (Fig. 10), but this process is well beyond the limits of applicability of LEFM, which does not take into account explicitly inelastic deformation. Note again that this fracture mode is possible only in the presence of stress concentrators, which can be related to the heterogeneities of the material (case analyzed in this work) typical for rocks or to structural heterogeneities also widely present in sedimentary basins. These heterogeneities and the related fracture initiation points are typical for natural joints and were extensively discussed in geological literature (e.g., McConaughy and Engelder, 2001; Pollard and Aydin, 1988). In Fig. 2c, the initiation point (zone) is

located at the place where the LVDT gage support was attached to the sample (see caption to this figure).

The numerical models thus predict a complication of the fracture process and roughening of the joint surfaces (walls) with increasing $P$. This is in total agreement with the experimental results reported in Chemenda et al. (2011a). If the stress level is still higher, the incipient shear bands (that do not evolve beyond the initiation stage during the described dilatancy banding) are resolved into macro shear bands and then shear fractures (Fig. 11). The angle $\psi$ between these fractures and $\sigma_1$ is rather high in this figure because the slope of the yield envelope adopted in the constitutive model changes abruptly at the BT-to-SD transition (Fig. 3a). In reality, the slope change is gradual
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Evolution of $\varepsilon^p$, of the octahedral strain rate $\dot{\gamma}$ (increment of the total octahedral strain during one calculation step), of the tensile strength $\sigma_t$, and of the mechanical state patterns in the model with Gaussian distribution of $\nu$ and a weak zone with $\sigma_t = 0.03$ MPa, $k_c = 0.3$ MPa. The initial stresses correspond to point 3 in Fig. 3a ($\sigma_2 = \sigma_1 = 0.6$ MPa). Only central segments of the images with $\varepsilon^p$ and $\dot{\gamma}$ patterns are shown, as beyond these segments, strains and strain rates are negligible.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 11. Evolution of the $\epsilon^p$ and mechanical state patterns in the model with Gaussian distribution of $\nu$ and a weak zone with $\sigma_t = 0.03$ MPa; $k_c = 0.3$ MPa. The initial stresses correspond to point 4 in Fig. 3a ($\sigma_2 = \sigma_1 = 0.66$ MPa).</figcaption></figure>

(Fig. 2a), which should be taken into account in future models. The impact of structural heterogeneities and notably those related to multilayer nature of the sedimentary basins (McConaughy and Engelder, 2001) should be addressed as well.

## Acknowledgements

This paper is invited in the frame of the prizes awarded by the 'Académie des sciences' in 2017. It has been reviewed and approved by Jean-Paul Poirier and by Vincent Courtillot

## References

* Bahat, D., 1991. Tectono-fractography. Springer-Verlag, Berlin, 354 p.

* Bésuelle, P., Desrues, J., Raynaud, S., 2000. Experimental characterisation of the localisation phenomenon inside a Vosges sandstone in a triaxial cell. Int. J. Rock. Mech. Min. Sci. 37 (8), 1223–1237.

* Chemenda, A.I., 2007. The formation of shear-band/fracture networks from a constitutive instability: theory and numerical experiment. J. Geophys. Res. 112, B11404, http://dx.doi.org/10.1029/2007JB005026.

* Chemenda, A.I., Mas, D., 2016. Dependence of rock properties on the Lode angle: experimental data, constitutive model, and bifurcation analysis. J. Mech. Phys. Solids. 96, 477–496.

* Chemenda, A.I., Nguyen, S.-H., Petit, J.-P., Ambre, J., 2011a. Mode I cracking versus dilatancy banding: experimental constraints on the mechanisms of extension fracturing. J. Geophys. Res. Solid Earth 116, no. B4.

* Chemenda, A.I., Nguyen, S.-H., Petit, J.-P., Ambre, J., 2011b. Experimental evidences of transition from mode I cracking to dilatancy banding.
C. R. Mecanique 339, 219–225, http://dx.doi.org/10.1016/j.crme.2011.01.002.

Fortin, J., Stanchits, S., Dresen, G., Guéguen, Y., 2006. Acoustic emission and velocities associated with the formation of compaction bands in sandstone. J. Geophys. Res. 111, B10203, http://dx.doi.org/10.1029/2005JB003854.

Ingraham, M.D., Issen, K.A., Holcomb, D.J., 2013. Response of Castlegate sandstone to true triaxial states of stress. J. Geophys. Res. 118, 536–552, http://dx.doi.org/10.1002/jgrb.50084.

Issen, K.A., Rudnicki, J.W., 2000. Conditions for compaction bands in porous rock. J. Geophys. Res. Solid Earth 105 (B9), 21529–21536.

Itasca, 2012. FLAC3D, Fast Lagrangian Analysis of Continua in 3 Dimensions, Version 5.0. Itasca Consult. Group, Minneapolis, MN, USA.

Lawn, B., 1993. Fracture of Brittle Solids, 2nd ed. Cambridge University Press, New York, 378 p.

Li, X., Feng, F., Li, D., Du, K., Ranjith, P.G., Rostami, J., 2018. Failure characteristics of granite influenced by sample height-to-width ratios and intermediate principal stress under true-triaxial unloading conditions. Rock Mech. Rock Eng. 51 (5), 1321–1345.

Ma, X., Rudnicki, J.W., Haimson, B.C., 2017. Failure characteristics of two porous sandstones subjected to true triaxial stresses: Applied through a novel loading path. J. Geophys. Res. Solid Earth 122, http://dx.doi.org/10.1002/2016JB013637.

Mas, D., Chemenda, A.I., 2015. An experimentally constrained constitutive model for geomaterials with simple friction–dilatancy relation in

brittle to ductile domains. Int. J. Rock. Mech. Min. Sci. 77, 257–264, http://dx.doi.org/10.1016/j.ijrmms.2015.04.013.

McConaughy, D.T., Engelder, T., 2001. Joint initiation in bedded clastic rocks. J. Struct. Geol. 23 (2–3), 203–221.

Nguyen, Si-H., Chemenda, A.I., Ambre, J., 2011. Influence of the loading conditions on the mechanical response of granular materials as constrained from experimental tests on synthetic rock analogue material. Int. J. Rock Mech. Min. Sci. 48, 103–115, http://dx.doi.org/10.1016/j.ijrmms.2010.09.010.

Pollard, D.D., Aydin, A., 1988. Progress in understanding jointing over the past century. Bull. Geol. Soc. Am. 100, 1181–1204.

Rudnicki, J.W., Rice, J.R., 1975. Conditions for the localization of deformation in pressure-sensitive dilatant materials. J. Mech. Phys. Solids 23, 371–394.

Savalli, L., Engelder, T., 2005. Mechanisms controlling rupture shape during subcritical growth of joints in layered rocks. Geol. Soc. Am. Bull. 117 (3–4), 436–449.

Vardoulakis, I., Sulem, J., 1995. Bifurcation analysis in geomechanics. Blackie Academic and Professional, New York, 462 p.

Wong, T.F., David, C., Zhu, W., 1997. The transition from brittle faulting to cataclastic flow in porous sandstones: Mechanical deformation. J. Geophys. Res. 102 (B2), 3009–3025.

Woodworth, J.B., 1896. The fracture system of joints, with remarks on certain great fractures. Boston Natl. Hist. Proc. 27, 163–183.