Determining crack initiation stress in unconventional shales based on strain energy evolution

Determining the crack initiation stress (Ci) for unconventional shale rocks is of critical importance in describing the entire failure process of unconventional shale reservoirs. We propose a new method to identify Ci values based on triaxial failure tests on four organic shale samples, attempting to improve the shortcomings of other methods. The new method is based on the relationship between crack development and strain energy evolution (SEE). Additionally, the proposed SEE method is compared with three widely used methods, including crack volumetric strain (CVS), moving point regression (MPR) and the lateral strain response (LSR), intending to examine the performance of different methods. The contrastive results indicate that the LSR method cannot determine Ci when the rock ruptures without volumetric dilatancy, which frequently occurs in the compression process of organic shales. Ci values obtained using the SEE method are consistent with those from the CVS and MPR methods. However, the proposed SEE method with a solid physical basis is more objective and stable than the CVS and MPR methods. The proposed method, from one aspect, compensates for the shortcomings of other methods when facing different failure modes in organic shales. From the other aspect, it provides a way to precisely determine Ci values for applications in wellbore stability evaluation and hydraulic fracturing design.


Introduction
Brittle failure and instability of unconventional formation rocks are critical for many geoengineering applications, such as underground excavation, tunnel spalling, borehole breakout and hydraulic fracturing (Martin et al. 1999;Diederichs 2007;Martin & Christiansson 2009;Damjanac & Fairhurst 2010;Zhao et al. 2016;Meng et al. 2020;Faraji et al. 2021). The progressive failure of rock materials is frequently induced by initiation, development and coalescence of preexisting or newborn cracks, which can be monitored using acoustic emission and microseismic analysis (Martin & threshold between the elastic deformation and crack initiation stages, is frequently used to characterise the onset of newborn cracks. This experimental study aims to examine a new method of determining Ci values for unconventional shale samples. In past decades, various methods have been proposed to detect Ci by stress-strain curves. Brace et al. (1966) have suggested that the onset of the nonlinear part in stress-strain curves is where new cracks are generated in the process of triaxial compression. Eberhardt et al. (1998) stated that Ci values of Lac du Bonnet granite could be determined from evolutions of the axial or volumetric stiffness with the applied triaxial stresses. However, these two methods are based on observations by the naked eye, which are too subjective to precisely select the deflected point in the measured curves. Martin & Chandler (1994) introduced a crack volumetric strain (CVS) method to determine Ci values for Lac du Bonnet granite. This method requires precise measurements of elastic mechanical parameters, such as Young's modulus and Poisson's ratio. Subsequently, Nicksiar & Martin (2012) updated an objective approach to obtain Ci values with lateral strain response (LSR) differences for low-porosity crystalline rocks. This method overcomes some flaws induced by subjective judgment or elastic parameters selection, but its hidden physical mechanism is still questionable (Zhang et al. 2021). Wen et al. (2018) have proposed a method based on the relative compression strain response to determine Ci for diverse rocks. The advantage of this method is avoiding users' judgment subjectivity.
However, the existing methods mentioned frequently focus on igneous or metamorphic rocks. The applicability of these methods to other rock types with more complicated structure, such as sedimentary rocks, is uncertain. Ning et al. (2018) applied these methods of identifying Ci to coal samples and compared with a method based on energy dissipation. Li et al. (2020) determined Ci values for anisotropic organic shales by calculating the volumetric strain of crack. Zhang et al. (2021) extended the Ci determination methods based on stress-strain curves to siliceous siltstone and limestone and obtained good results after comparing with the acoustic emission results.
Although many tentative methods have been proposed to determine where new cracks are initiated for various rocks, no individual method has been proved to be the best way to determine Ci values. Additionally, most of the methods are too subjective and short of a solid physical basis. In contrast to crystalline rocks, the stress-strain curves of unconventional shales become more complicated. It is challenging to precisely determine Ci values by solely relying on the stress-strain curve morphology.
The nature of triaxial fracturing is the energy conversion. In this work, we aim to objectively identify Ci values for unconventional shale samples based on SEEs. The perfor-mance of the proposed method is compared with that of three widely used methods. First, we introduce the selected shale samples and the apparatus used for triaxial compression failure tests. Second, the stress-strain behaviors of four shale samples are analysed. Then, four methods of determining Ci values are displayed and discussed. Based on the analysis, we emphasise the reliability of the new method in future applications.

Sample description
Four organic cores with obvious bedding planes are extracted from a shale oil reservoir located in north China. The original burial depth of four full-diameter cores ranges from 2000 to 2600 m. Subsequently, four cylindrical samples are drilled from the full-diameter cores in the direction perpendicular to bedding planes. Four cylindrical samples are named F-2, F-4, Z-11 and Z-19, as shown in figure 1. The ends of each plug are machined to be parallel to the diameter direction with a tolerance of ±2.5 × 10 -3 cm. The length-diameter ratio is about 2:1 based on the International Society for Rock Mechanics (ISRM's) suggestion for rock mechanics experiments. Then, cylindrical samples are placed in an electric oven for drying with a constant temperature of 80ºC until the weight does not change. The bulk density for each dry rock is calculated from the rock weight and volume. A helium porosimeter is used to measure the grain density. Subsequently, the porosity of each sample is derived based on the bulk and grain density, as shown in Table 1.
Additionally, we carry out X-ray diffraction analysis on F-4 and Z-19 samples to explore the mineral compositions. F-4 shale is mainly composed of 36.2% quartz, 15.9% feldspar and 47.8% clay. Z-19 shale is primarily made up of 41.2% quartz, 31.1% feldspar and 15.4% clay. In addition, the total organic carbon in the F-4 shale (1.6%) is higher than that in the Z-19 shale (0.4%). From the thin section image in figure 2, taking from the end face of F-4 shale, solid minerals (e.g. quartz, feldspar) are discretely distributed in soft clay minerals. Black kerogen is scattered across clay minerals.

Experimental system and procedure
We conduct triaxial compression failure tests on four shale samples using the MR-RM 3000 rock mechanics testing system. All the measurements are carried out at room-dry conditions. The axial load is applied using a ball screw press, which uses a microprocessor-based electronic servo drive to provide high precision load and position control. This system offers an internal load of up to 60 000 pounds and a confining pressure of ∼70 MPa. The sample assembly inside 643 Figure 1. Pictures of ruptured samples after triaxial compression tests.  the confining vessel is illustrated in figure 3. The cylindrical sample is covered with a thin heat-shrink Teflon tube and placed between the top and bottom platens. There is an internal load cell under the bottom platen to record the deviatoric stress along with the measurements. The axial strain is measured with a couple of linear variable differential transformers fixed between the top and bottom platens. The lateral displacement is measured with two pairs of strain gauges attached to a cantilever bridge. The average of lateral strains in two orthogonal directions is used for the following analy- sis. The volumetric strain is defined as axial strain plus twice the lateral strain . More details of the experimental apparatus can be found in Wang et al. (2020a). After closing the confining vessel, the confining pressure is raised to ∼20 MPa with a rate of 0.5 MPa s −1 . The 20 MPa confining pressure, to some extent, simulates the 644 Figure 4. Stress-strain curves for four measured shales.
in situ effective horizontal stress according to the samples' burial depths. Subsequently, the deviatoric stress is applied when the piston moves downward to touch the top platen. The deviatoric load compresses the rock sample along the axial direction at a constant confining pressure of ∼20 MPa until it ruptures. During the triaxial compression measurements, the axial strain rate is controlled at 3 × 10 -6 s. Figure 4 displays the stress-strain curves of four shale samples. Overall, all strains respond to the applied deviatoric stress in nonlinear patterns. Z-19 shale exhibits the highest compressive strength compared to other shales (F-2, F-4, and Z-11). F-2 and Z-19 shales display apparent volumetric strain reversals before rock failure, whereas it is difficult to determine volumetric strain reversals for F-4 and Z-11 shales. The evolutions of strains with the applied compression, to some extent, describe the processes of crack nucleation, crack growth and crack interaction, which results in macroscopic fracture under compression (Amann et al. 2011). As stated by Brace et al. (1966), the micro-crack development processes can be depicted in five stages: pre-existing micro-crack closure, linearly elastic deformation, new crack initiation, crack damage and failure. The corresponding four key stress thresholds are pre-existing micro-crack closure stress (Cc), new crack initiation stress (Ci), crack damage stress (Cd) and compressive failure strength (Cp) (Brace et al. 1966;Bieniawski 1967;Martin & Chandler 1994;Li et al. 2020). Crack damage stress represents the threshold transiting from a stable crack growth to unstable crack damage, which can be acquired from the reversal of the volumetric stress-strain curve. From this viewpoint, it is easy to determine the crack damage stress for F-2 and Z-19 shales but difficult for F-4 and Z-11 shales. In addition, there is a linearly elastic deformation stage in stress-strain or strain-strain curves, which is used to derive static Young's modulus and Poisson's ratio according to the ISRM standard (Wang et al. 2020b), as shown in Table 1. By comparison, Z-19 shale with the lowest porosity has the largest Young's modulus and Poisson's ratio.

Stress-strain evolutions
One challenge in investigating crack development processes during the compression failure test is to precisely identify the crack initiation (Li et al. 2020;Zhang et al. 2021). Many researchers have proposed a variety of methods, including CVS, Martin & Chandler 1994), the linear strain response (LSR, Nicksiar & Martin 2012), moving point regression (MPR, Eberhardt et al. 1998) and acoustic emission (AE, Zhang et al. 2021), attempting to overcome these challenges. This study proposes a new approach based on strain energy evolutions (SEE) by analysing stress-strain curves, aiming to identify accurate Ci. Meanwhile, three widely used methods (CVS, LSR and MPR) are also applied to detect Ci 645 Figure 5. The strain energy density at point X on the stress-strain curve of Z-11. using stress-strain data. The reliability of each method is investigated by comparing Ci results from different methods.

Proposed strain energy evolution (SEE) method
In triaxial compression tests, with the deformation of rocks, energy input, accumulation and release occur concurrently. At the initial stage of compression, the rock specimen absorbs and stores a mass of elastic strain energy charged by the external load. Thereafter, a part of the energy is dissipated to initialize and develop cracks. With the increase of the dissipated strain energy, the structural collapse is induced by more and more damage. Hence, the crack development is strongly related to the dissipated strain energy. Ning et al. (2018) proposed that the energy dissipation derived from the cyclic loading test could describe crack growth. However, the cyclic loading test is expensive and time-consuming due to its complex procedures. Moreover, the cyclic loading with different stress magnitudes will alter the original rock structure, which subsequently influences rock mechanical properties. Hence, we calculate the SEE along with monotonous compression failure tests based on a simplified model, as shown throughout equations (1-4).
In the monotonous compression test, the derived static Young's modulus in Table 1 serves as an input of the model to obtain the elastic and dissipated strain energy density (U e and U d ) (Xie et al. 2009;Tarasov & Potvin 2013;Munoz et al. 2016). The total strain energy density (U t ) is calculated as the area below the curve between the deviatoric stress and the axial strain (figure 5), which represents the rock's energy status. Therefore, we calculate different kinds of strain energy at any specific stress level (figure 5). According to the first thermodynamics law (Clausius 1857), strain energy conversions in rocks can be expressed as where E represents the static Young's modulus with a unit of GPa. i is the deviatoric stress at point i of the stress-strain curve with a unit of MPa, and i is the corresponding axial strain in response to i . According to the equations, we calculate the SEE as a function of the deviatoric stress for four shale samples, as shown in figure 6. During stages of the crack closure and elastic deformation, no new cracks are created and developed. The curve of the elastic energy density almost overlaps with that of the total strain energy density, indicating that nearly no energy is dissipated during these two stages. Subsequently, the density evolutions between the elastic energy and the total strain energy begin to divaricate when the deviatoric stress reaches a certain level, which indicates the occurrence of the dissipated strain energy. The dissipated strain energy contributes to the new crack growth, suggesting that the rock enters the nonlinear stage. Thus, the deviatoric stress, where the forked point appears (figure 6), is considered Ci. Ci coincides with where the dissipated strain energy occurs. By contrast, in figure 6, Z-19 shale with the lowest porosity exhibits the largest Ci value.

Crack volumetric strain (CVS) method
The CVS has proved to be an effective way to determine the crack initiation stress, as stated by Martin & Chandler (1994) after performing a suite of compression failure tests on various granite samples. The CVS (V c ) is defined as the difference between the measured volumetric strain (V) and the elastic volumetric strain (V e ). It can be expressed as where E and are the static Young's modulus and Poisson's ratio, respectively, d is the deviatoric stress with a unit of MPa, and a and l are the axial and radial strain, respectively. Figure 7 exhibits the calculated V c as a function of the axial strain for four shale samples. Initially, the CVS almost remains constant values approaching zero, except for Z-11 shale. Once a threshold axial strain is reached, the CVS 646 Figure 6. The calculated strain energy density (i.e. total, elastic and dissipated) as a function of deviatoric stress for four shale samples. decreases (F-2 and Z-19) or increases (F-4 and Z-11) with the continued compression. This threshold axial strain corresponds to where cracks are initially generated, namely Ci. The different trends when the deviatoric stress is beyond Ci are attributed to the magnitude relations between axial and radial strains. The continuous increase of the CVS for F-4 and Z-11 is caused by the lack of dilatancy of the volumetric strain curve, as shown in figure 4, which denotes larger axial deformation, in contrast, to double of the radial deformation. However, regardless of the increasing or decreasing trend, the determination of Ci is not affected. The relatively flat variation of the CVS before Ci indicates that the elastic deformation mainly contributes to the total volumetric strain. Thus, as thresholds between elastic and plastic stages, Ci values can be detected in figure 7.

Moving point regression method (MPR)
Another effective way to identify Ci values is through determinations of the volumetric stiffness (Eberhardt et al. 1998). With the MPR method, the slope of the volumetric stressstrain curve in a selected window is derived, which is defined as the volumetric stiffness. Figure 8 exhibits the evolutions of the volumetric stiffness with the deviatoric stress for four samples. Generally, the volumetric stiffness evolves with an increasing or decreasing trend before reaching the crack closure stress threshold. With continued load, the stiffness levels off, indicating that the linearly elastic stage is reached. The mutation of the volumetric stiffness after the elastic stage indicates the crack initiation. For F-2 and Z-19, the volumetric stiffness keeps increasing until the reversal of the volumetric stress-strain curve is reached. However, there is no volumetric dilatancy for F-4 and Z-11 shales during the failure tests, resulting in a continuously decreasing trend of the volumetric stiffness. When the deviatoric stress reaches the failure strength, the stiffness abruptly shifts from positive to negative, indicating that the rock is ruptured. From figure 8, Ci values for all the measured samples are determined when the volumetric stiffness starts to change from a relatively consistent value in the elastic stage.

Lateral strain response method (LSR)
LSR method is another way to detect Ci from stress-strain curves (Nicksiar & Martin 2012). The lateral strain evolution with the deviatoric stress is examined from zero stress to the crack damage stress. Meanwhile, a reference line is selected from zero stress to the crack damage stress, as shown in figure 9. Thus, the difference between the reference line and the lateral strain curve at the same deviatoric stress is 647 Figure 7. The calculated crack volumetric strain as a function of axial strain for four shale samples. defined as the lateral strain difference (figure 9). Subsequently, the derived lateral strain differences are fitted polynomially. The maximum strain difference on the best fitting parabola indicates Ci. It is pertinent to mention that there are no volumetric reversal points for F-4 and Z-11 shale, as denoted in figure 4. Hence, we use the peak strength instead of the crack damage stress for calculating the lateral strain differences. All the results from the LSR method are displayed in figure 10.

Ci determination for two groups of shales
Based on stress-strain curves, Ci values are determined for all the measured samples using four methods, as shown in Table 2. Meanwhile, the mean value and standard deviation for each sample are also presented. Overall, Z-19 shale with the smallest porosity shows the largest Ci regardless of the used method. Four samples can be categorised into two groups based on the morphology of volumetric stress-strain curves in figure 4. The first-group samples (F-2 and Z-19) experience obvious volumetric dilatancy once the sample approaches failure. For F-2 shale, Ci values determined using four methods are similar with an average of 47.4 MPa. The SEE and LSR methods display Ci values as 77.4 and 74.8 MPa for Z-19 shale, respectively, which are higher than those based on the CVS and MPR methods. The second group, including F-4 and Z-11, represents the shale rupture without volumetric dilatancy. For F-4 and Z-11, Ci values range from 39.7 to 70.5 MPa. The LSR method estimates an anomalously larger Ci value than other methods. There is no noticeable difference for the Ci value obtained from the SEE, CVS and MPR methods.
For rocks without volumetric dilatancy, the rock volume is compacted accompanied by ductile damage under compression (Scott & Nielsen 1991;Popp et al. 2001). Meanwhile, based on the shale failure model, in the vertical shales with the loading direction perpendicular to bedding planes, soft layers lead to significant displacement in the axial direction but rare lateral strain change (Wang et al. 2020a). Hence, the lateral strain is extremely smaller than the axial strain, giving rise to the continuous increase of the volumetric strain in the positive direction. However, the LSR method, based on the lateral strain response, will not be sensitive enough when the lateral strain is too small, as exampled by F-4 and Z-11 samples in figure 9. Therefore, the anomaly Ci is observed for F-4 and Z-11 from the LSR method.
Additionally, the unconventional reservoir shales generally contain bedding layers and their structures are characterised by the intrinsic anisotropy (Wang et al. 2021). The underground shales are also subject to complex in situ stress environments. The bedding layer orientation with respect to the in situ stress would significantly influence the mechanical properties of shales (Li et al. 2020;Wang et al. 2020a). However, as demonstrated by Li et al. (2020) after fracturing 648 Figure 8. Evolutions of the volumetric stiffness with applied deviatoric stress for four shale samples. Figure 9. Lateral strain as a function of deviatoric stress ranging from zero to crack damage stress for F-2. The method to calculate the lateral strain difference is displayed.
Longmaxi shales with seven different bedding orientations, the determined Ci values are independent on the bedding orientation. They attributed the crack initiation to the randomly distributed microcracks due to the complex geological processes. But the bedding orientation would greatly influence the crack propagation. The initiation of microcracks in the process of triaxial compression is the focus of this research. Therefore, the presence of beddings in the vertical shales would not influence the determination of Ci.

Evaluation of different methods for Ci determination
For the CVS method, it is relatively effective to establish Ci for most of the measured samples. But the Ci value from this method is lower than that from other methods for Z-19 shale. The Ci value is strongly related to the precision of elastic properties according to equations (1)-(4). Usually, the measurement of Young's modulus is trustworthy in the compression test. However, it is difficult to obtain an accurate Poisson's ratio, especially for shales with pronounced anisotropy (Wang et al. 2021). Therefore, the result is significantly affected by Poisson's ratio in this method (Eberhardt et al. 1998;Zhang et al. 2021).
For the MPR method, consistent performance is observed in all the measured samples. Compared with other methods, this method obtains slightly lower Ci values for each tested sample. Nevertheless, this method is too subjective to discover the sudden change of volumetric stiffness after the elastic stage. The measured stress-strain data and the selected window length induce fluctuations of the volumetric stiffness curve (figure 8), which enhances the difficulty of determining Ci precisely.
The LSR method is reliable for determining Ci for F-2 and Z-19, which exhibit volumetric reversal as the rock approaches failure. However, it fails to determine the crack damage stress once no volumetric dilatancy occurs, such as F-4 and Z-11 shales. The LSR method is strongly influenced by the determination of the crack damage stress, although it 649 Figure 10. The lateral difference as a function of deviatoric stress for all measured samples. is less subjective than the MPR method (Zhao et al. 2015). Additionally, the physical meaning of this method is still ambiguous (Zhang et al. 2021). Hence, the LSR method is not suitable to identify Ci for the shale, which does not display volumetric dilatancy during the failure process. The proposed SEE method displays a robust performance to detect the Ci values based on the consistent results for each of the shale samples. This method can be used not only for the conventional failure mode (F-2 and Z-19) but also for the non-dilatancy failure mode (F-4 and Z-11). Compared with the MPR method, it is more objective at detecting Ci.
Meanwhile, different from the CVS method, it only requires Young's modulus and the stress-strain data as inputs, which avoids the potential error induced by Poisson's ratio. It should be noted that the SEE after elastic damage might not represent the precise evolution process due to the lack of the cyclic loading test. However, it does not affect the determination of the threshold between the elastic and plastic stages if accurate Young's modulus is provided.
Precisely determining Ci values is of great significance in evaluating the progressive failure of rocks. The Ci value could be used for determining the necessary parameters in establishing the Hoek-Brown failure criterion (Martin et al. 1999;Andersson et al. 2009;Wang et al. 2020b). In addition, the distance between the crack initiation and the rock failure could be used as a brittleness index, which further guides the optimal design of hydraulic fracturing for unconventional shale reservoir development (Stanchits et al. 2014;Wen et al. 2018;Chen et al. 2019). The pre-existing methods are mostly limited by the human subjective judgment or the accuracy of elastic parameters and cannot determine where the microcracks are initiated precisely. By comparing, the proposed SEE method provides an effective and stable way to identify Ci for unconventional shales on the basis of energy conversion. 650

Conclusion
Methods to detect the crack initiation stress (Ci) are discussed through triaxial compression failure tests on four organic shale samples. Due to the close link between strain energy and crack development, the SEE method is developed to identify Ci in this study. Based on the results from different methods, the CVS method and the MPR method display some shortcomings, which can induce errors in the Ci determination. In addition, the lateral strain method (LSR) overestimates the Ci value for the sample without volumetric dilatancy under the failure process. The proposed SEE method is considered stable and objective to obtain accurate Ci for all the tested samples because it avoids the shortcomings of other methods. It is appliable in different failure modes and practical for Ci determination based on Young's modulus and stress-strain data. Therefore, the SEE method provides a new viewpoint to determine Ci accurately for engineering applications in shale reservoirs.