Effect of heterogeneity on failure of natural rock samples

A carbonate sample extracted from the depth of about 10 kft was subjected to uniaxial loading while the confining stress remained constant. Post-experiment inspection of the sample showed an inclined crack at an angle less than 20° to the horizontal. This subhorizontal crack orientation was contrary to the expected 45° inclination, the plane of the maximum shear stress. Coincidentally, as shown by CT-scan prior to loading, there was a boundary between two layers of different density inside the sample located almost exactly where the crack appeared. This density difference has arguably translated into the contrast in the elastic properties at the boundary. The hypothesis is that because of this elastic heterogeneity, an incipient crack developed at the boundary due to the unavoidable tensile stressing of the sample as it was brought to the benchtop from its original state of high confining stress at depth. Controlled uniaxial compression made the sample slip along this crack, which then developed into a prominent feature. This assumption was corroborated by a numerical experiment showing a strong von Mises stress concentration at the elastic contrast boundary during hydrostatic tensile loading. Another sample, from the same formation, but without strong density heterogeneity, exhibited a classic 45° crack after uniaxial loading. These results provide a novel and important insight into the mechanics, breakage, and strength of natural rock.

Tests and results. Carbonate samples C23 and C8 extracted from a well drilled in a Saudi oil field were cleaned and dried in the laboratory. Their dynamic Poisson's ratio ν and Young's modulus E were computed from the bulk density ρ b and the P-and S-wave velocities, V p and V s , measured at 30 MPa hydrostatic stress as and are listed, together with the depth of their original location, in Table 1. www.nature.com/scientificreports/ Both samples were first hydrostatically loaded to 30 MPa. Afterwards, the vertical stress was gradually increased until failure, with the confining stress remaining constant. Finally, after a period of post-failure deformation, the samples were axially unloaded to the initial hydrostatic stress 30 MPa. All tests were performed using a load-frame Autolab-1500 supplied by New England Research. The elastic-wave velocities were measured on room-dry samples using ultrasonic transducers with 750 kHz central frequency. The axial and radial deformations were measured by linear variable differential transformers (LVDT). The results of these tests are shown in Fig. 2.
In both tests, observed was a textbook behavior with almost linearly elastic deformation followed by failure, plastic deformation, and elastic unloading. The unloading curves essentially parallel the loading curves. The post-yield flow time for C23 was 30 min, while it was 8 min for C8. This explains why the post-yield deformation in C23 was about three times that in C8.
First sample: post-test and pre-test images. Post-test CT-scan images (vertical slices of the original 3D images taken along the central axis) for sample C23 are shown in Fig. 3. The major crack in the lower half of the sample is inclined at about 20° to the horizontal and surrounded by minor cracks at the circumference of the sample.

Methods
The hypothesis is that an incipient subtle crack developed at the boundary as the sample was brought from its natural depth of about 10 kft to the ambient conditions at the benchtop. This crack became the "plane of weakness" along which the sample failed during the experiment.
Indeed, the original confining stress at depth was on the order of 30 to 40 MPa. Upon its extraction, the sample was inadvertently subjected to tensile hydrostatic stress as the stress was reduced from its original value to zero. As a result, an incipient crack appeared at the boundary inside the sample and then developed into a major fracture during the uniaxial loading.
To test this hypothesis, a digital cylindrical sample with a plane boundary cutting it at about 20° to the horizontal was created. The Young's moduli E above and below this discontinuity were assigned 56 and 36 GPa, respectively, while the Poisson's ratio ν was the same 0.3.
The density distribution inside the model was assigned based on the CT-value contrast shown in Fig. 3. The total porosity was computed from density and using the mass-balance law by assuming that the mineral is pure calcite with density 2.71 g/cc. Finally, this porosity contrast was translated into the elastic-moduli contrast by assuming the stiff-rock effective-medium model described in, e.g., 10 . This model accurately describes the velocity-porosity wireline data from the well where the samples were extracted from (see 11 ).
Uniform hydrostatic tension of 30 MPa was applied to this digital object and the von Mises stress was computed using the finite elements method (FEM). These computations were done using a commercial FEM package COMSOL, where the digital model was meshed by triangular elements fine enough to delineate the boundary of the elastic contrast and resolve the ensuing stress concentration at this boundary. The results shown in Fig. 4 indicate a sharp σ VM concentration at the boundary between the two domains with varying E . This is where the initial crack arguably developed as the sample was lifted from its original depth. Notice also the von Mises stress concentration at the circumference of the elastic discontinuity. Their location is where the minor cracks surrounding the main crack developed.
The stresses shown in Fig. 4 are in MPa, although any units could be used with the respective stress boundary conditions since the simulations were conducted under a linear elasticity assumption.
Of course, the same stresses but with the opposite sign would develop under hydrostatic compression. However, because rocks are weaker in tension than they are in compression (e.g., 12 ), still valid is the hypothesis that the "plane of weakness" in this sample developed during its extraction from depth, prior to laboratory experiments.
Second sample: post-test and pre-test images. The post-and pre-test images for sample C8 are shown in Fig. 5. It also fractured, but unlike in the first sample, the crack is inclined at the textbook 45° to the horizontal.
There are also high-porosity darker strips visible in the pre-test image. However, unlike in C23, they are oriented at approximately 45° to the horizontal. Once again, an incipient crack might have developed along one of these discontinuities during sample recovery. Numerical tests, same as shown in Fig. 4, indicate that no www.nature.com/scientificreports/ matter how an elastic property discontinuity is oriented, sharp stress concentration appears along its trajectory under hydrostatic tension.
Third sample: post-test and pre-test images. The hypothesis put forward here is further supported by the testing and respective pre-and post CT-scan images of the third carbonate sample. During the test, this sample's behavior was qualitatively the same as that of the first two samples as shown in Fig. 2. A 3D CT-scan image taken after the test exhibited massive fracturing as shown in a 2D axial slice of the image in Fig. 6 (right). The locations of these fractures coincided with CT-intensity heterogeneity apparent in the pre-test image (Fig. 6, left). As in the first two samples, this intensity heterogeneity translated into density and elastic moduli heterogeneity, hence resulting in von Mises stress concentrations during tensile unloading of the core as it was evacuated from its original depth.

Discussion and Conclusion
The effect of structural heterogeneity on crack initiation was discussed in 13 , where numerical modeling showed that the micro-heterogeneity played an important role in controlling both the micromechanical behavior and the macroscopic response when subjected to uniaxial compression loading. The crack-initiation stress was found to be controlled primarily by the micro-scale geometric heterogeneity. Earlier experimental results (e.g. 14 ) also showed that microstructural heterogeneity plays a key role in creating local stress concentrations. Recently 15 ,  www.nature.com/scientificreports/ showed that grain boundary sliding can induce strain and rotation at the respective interfaces. In the new work reported here, a sample failure is observed as well, arguably due to tension-induced stress concentration along a structural heterogeneity. It is evident from the CT-scan images in Figs. 3, 5, and 6 that the samples under examination contain relatively large vugs. A question arises whether the failure of a sample could have occurred due to the stress concentration around these features. A visual inspection of the aforementioned images indicates that this is not the case. Indeed, the fracturing occurred at the density contrast boundaries rather that at the vugs. To quantitatively assess the possibility of fracturing due to stress concentration at the vugs, we present here a simple model of an elastic body with isolated, as well as interacting spherical and elliptical vug-like inclusions. Figure 7 shows that stress concentrations under tension are localized around isolated inclusions as would be also evident from classic closed-form solutions 16 . Where the inclusions interact, once again, stress concentrations are localized around and between the inclusions and are unlikely to produce a plane of weakness cutting through the entire sample.
This work followed the classic scientific method: from the observation to a hypothesis to its validation to additional testing. Specifically, in trying to understand reasons for the occurrence of a subhorizontal (rather than inclined at 45°) crack in a rock sample subject to uniaxial loading, it was noticed that the sample exhibited a clearly visible density and structural discontinuity located exactly where the crack appeared. This discontinuity most likely translated into an elastic moduli contrast. The hypothesis is that an incipient fracture was generated along this discontinuity as the sample was lifted from the depth and unavoidably unloaded during this extraction. Later, during loading, the sample failed along this subtle plane of weakness. A numerical experiment confirmed