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 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.

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

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 $\overline{\tau }$ is von Mises's stress; $R={\overline{\tau }}_c/{\overline{\tau }}_e$; ${\overline{\tau }}_c$ and ${\overline{\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

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

where μ_c, k_c, μ_β 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 μ_c and k_c as well as of the tensile strength σ_t (${\mu }_{\mathrm{c}}^0$, $k_{\mathrm{c}}^0$ and ${\sigma }_{\mathrm{t}}^0$) are defined from the GRAM experimental data (from the simplified σ₃(σ₁) curve in Fig. 3a) to be

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

where ${\sigma }_{\mathrm{t}}^0$ = 0.07 MPa; ${\varepsilon }_{\mathrm{t}}^0$ = 0.5 × 10⁻³ and ɛ⁰ = 1 × 10⁻² (${\varepsilon }_t^0$ <<ɛ⁰); ${\mu }_{\mathrm{c}}^{\mathrm{end}}$ = 1.7, ${\epsilon }_{\mathrm{t}}=\int \mathrm{d}{\epsilon }_3^{\mathrm{p}}$ and ${\epsilon }_{\mathrm{t}}={\epsilon }_{\mathrm{t}}+\int \mathrm{d}{\epsilon }_1^{\mathrm{p}}$ (the superscript “p” stands for plastic or inelastic). Note that σ₁ 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 β_c evolves (reduces) with ɛ in the same way as μ_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 ν = 0.25 and Young's modulus E = 6.7 × 10⁸ 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.

The modeling setup (Fig. 3b) is two-dimensional in the sense that the model thickness in the σ₂-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 σ₁ = σ₂ and to a sufficiently large σ₃ for the stress state to be above the initial yield envelope in Fig. 3a (this corresponds to the purely elastic mechanical state). Then, σ₃ is progressively decreased by applying velocities V₃ to the σ₃-normal model boundaries, while maintaining the initial values of σ₁ and σ₂ by applying velocities V₁ and V₂ 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₂ is not zero, the models are not plane-strain.

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

4.1 Point 1 in Fig. 3a, σ₁ = σ₂ = 0.3 MPa

4.1.2 “Homogeneously heterogeneous” model with small variations of the Poisson's ratio ν

The Gauss distribution of ν 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 σ₃ 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 σ₃. 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, SD mechanism being operated only within the weak zone (Figs. 6c, d).

4.2 Point 2 in Fig. 3a, σ₁ = σ₂ = 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 propagates due to both BT and SD mechanisms (Fig. 8b–d). The oblique fracturing process is followed by an almost σ₃-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 σ₃-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, σ₁ = σ₂ = 0.6 MPa

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, σ₃ is tensile. Its maximum magnitude reaches around 0.06 MPa, which is somewhat less than ${\sigma }_{\mathrm{t}}^0$ = 0.07 MPa. The nominal stress (which would be measured at the model (sample) end in the experimental test) is σ₃ = −0.0011 MPa, which is much less in magnitude than ${\sigma }_{\mathrm{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{\gamma }$ and σ_t patterns in Fig. 10. During generally σ₃-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 σ₃. In the next model, it is the oblique direction that is definitively followed by the deformation band.

4.4 Point 4 in Fig. 3b, σ₁ = σ₂ = 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 σ₃-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).

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 fractures, are supposed to form exclusively when the tensile stress σ₃ reaches the tensile strength σ_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 }_{\mathrm{t}}^0$), joints can form at |σ₃| < <${\sigma }_{\mathrm{t}}^0$ (Figs. 8 and 9). Here σ₃ 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 σ₃, which can be measured experimentally in the laboratory. At the micro-scale, |σ₃| is tensile, but lower than ${\sigma }_{\mathrm{t}}^0$. Purely tensile, σ₃-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 σ₃-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 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 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 ψ between these fractures and σ₁ 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 (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.