An Evaluation Index of Fracability for Reservoir Rocks Based on Fracture Process Zone

A reliable evaluation method for the fracability (i.e., ability to generate abundant cracks) of reservoir rocks is a critical issue for maximum hydraulic fracturing efficiency. Most previous fracability indices lacked enough rationality and practicability, and thus could not consistently provide a reliable evaluation. We suggest a new fracability index called crack tolerance, which is represented by the maximum radius of the fracture process zone at the crack tip of a cracked chevron notched Brazilian disk specimen, which corresponds to the critical state for unstable propagation of the notched crack. In experiments and simulations based on the discrete element method, we showed quantitative methods to conveniently determine the value of the crack tolerance and showed that specimens with a greater crack tolerance generated more cracks before rupture and had complex morphologies, which would indicate stronger fracability. The crack tolerance can well characterize the effects of structural and loading conditions, including the grain size heterogeneity, bedding orientation, and environmental temperature, on fracability, and the inherent heterogeneity of rock is the physical basis for it as a fracability evaluation index. Our studies showed the rationality and practicability of this index and provide hints for how to produce abundant complex cracks in reservoirs.


Introduction
Human beings' mining engineering and energy resource exploitation extensively involve the generation and propagation of cracks within rock materials. Hydraulic fracturing is widely used to enhance the fluid conductivity of reservoirs of oil, gas, and geothermal resources. A reliable evaluation of the rock fracability (i.e., ability to generate abundant cracks) is important for hydraulic fracturing [1].
Brittleness, which is generally viewed as a property (or ability) of solid material that ruptures with little appreciable permanent deformation, has long been considered approximately equivalent to fracability, because it shows empirical relevance to the possibility of crack propagation: reservoir comprising brittle rocks usually responds well to stimulation, whereas preexisting and hydraulic fractures tend to heal rather than to propagate in a less brittle reservoir. This is probably attributed to less energy consumed by the ductile deformation of brittle rock materials [2].
In the past decades, a variety of brittleness indices have been developed to evaluate its effect [3,4], which can be classified into several broad categories: (1) Based on mineral composition (e.g., [5]), especially the weight or volume proportion of hard minerals such as quartz: a positive correlation seems to exist between the brittleness and mineral contents of rocks. However, such indices do not consider many other factors that also contribute to brittleness, such as grain size and loading conditions. (2) Based on elastic parameters (e.g., [6]): for example, rocks with a large Young's modulus and small Poisson's ratio are assumed more brittle. However, such indices can be controversial because many laboratories and in situ observations [7,8] contradict this assumption. (3) Based on strength: for example, one such index is the ratio of tensile and compressive strengths [9]. Such indices are easily measured, but they lack a physical correlation to brittleness and cracking propagation mechanisms. Thus, these indices may return similar values for various types of rocks with different levels of brittleness. (4) Based on characteristics of the stress-strain curve such as the relative stress drop, post-peak modulus, and various combinations [10,11]: these indices characterize rock brittleness well and are widely used in predicting the rockburst proneness. However, high brittleness does not consistently represent strong fracability because brittle rock can also act as a barrier to hydraulic fracturing [12].
In summary, many brittleness indices currently popular in fracability evaluation for reservoir lacks mechanical relevance to the rock cracking process. On the other hand, the evaluation indices used in other areas (e.g., those used to estimate rock cuttability [13]) are usually inapplicable for reservoir fracability evaluation owing to the essential differences of physical meaning between brittleness and fracability. Thus far, few evaluation indices of rock fracability meet the following requirements [3]: (1) Has a firm physical basis; (2) Consider the heterogeneity of rock material; (3) Be convenient to measure; (4) Characterize the effects of loading conditions.
To address this issue, we propose a new evaluation index for rock fracability that we call the crack tolerance. See Section 2 for its definition. Sections 3 and 4 show the experimental measurement of this new index and the corresponding numerical simulation results, respectively, to demonstrate the rationality of the index. Based on these analyses, the effects of several characteristics of the rock materials on the crack tolerance are discussed in Section 5. This study demonstrated the physical rationality of the crack tolerance as an evaluation index and analyzed the effects of the rock structure and loading conditions on the crack tolerance in an effort to extend our understanding of rock fracability and provide hints for how to produce more cracks in the reservoir.

Fracture Process Zone and Crack Tolerance
Numerous researchers have revealed that the propagation of macroscopic cracks within rock under tension is attributed to progressive generation, interaction, and nucleation of micro-cracks from the macroscopic crack tips as follows [14]. When the imposed tensile load is small, only a few independent micro-cracks can arise around each crack tip ( Figure 1a). As the tensile load increases, the distribution range of the micro-cracks expands, and their density increases. They interact with each other and coalesce (Figure 1b) to cause a gradual macroscopic propagation of the preexisting crack (Figure 1c,d). These micro-cracks indicate nonlinear deformation in the region around a crack tip preceding crack unstable propagation, which is referred to as a fracture process zone (FPZ) [14]. Crack propagation in tensile mode is most common in hydraulic fracturing because the effect of hydraulic pressure imposed on the crack surface approximates remote tensile stress in nature; additionally, rocks have a much lower tensile strength than compressive and shear strengths. Thus, cracks easily propagate driven by an injected fluid. The principal stresses at a tensile crack tip can be described as [15] where σ1 and σ2 are maximum and intermediate principal stresses, KI is tensile stress intensity factor, r and θ are polar radius and polar angle for polar coordinate system from the tip. Note that the minimum principal stress not listed here equals to zero. The range of FPZ (i.e., its size) is calculated based on the hypothesis that nonlinear deformation occurs within a region around crack tip when the local stress state satisfies a certain criterion (e.g., tensile strength criterion for rock materials, von Mises criterion for metal materials). The tensile cracks are assumed to propagate parallel to their own plane (i.e., θ = 0) when the σ1 reaches the tensile strength of the rock (σt), because the critical state of crack propagation is attained, which corresponds to the maximum size of the FPZ: t = IC √2π c cos 0 (1 + sin 0) which leads to Crack propagation in tensile mode is most common in hydraulic fracturing because the effect of hydraulic pressure imposed on the crack surface approximates remote tensile stress in nature; additionally, rocks have a much lower tensile strength than compressive and shear strengths. Thus, cracks easily propagate driven by an injected fluid. The principal stresses at a tensile crack tip can be described as [15] where σ 1 and σ 2 are maximum and intermediate principal stresses, K I is tensile stress intensity factor, r and θ are polar radius and polar angle for polar coordinate system from the tip. Note that the minimum principal stress not listed here equals to zero. The range of FPZ (i.e., its size) is calculated based on the hypothesis that nonlinear deformation occurs within a region around crack tip when the local stress state satisfies a certain criterion (e.g., tensile strength criterion for rock materials, von Mises criterion for metal materials). The tensile cracks are assumed to propagate parallel to their own plane (i.e., θ = 0) when the σ 1 reaches the tensile strength of the rock (σ t ), because the critical state of crack propagation is attained, which corresponds to the maximum size of the FPZ: which leads to r c = (K IC /σ t ) 2 /(2π),

of 18
where K IC is the tensile fracture toughness, and r c is the maximum FPZ size. In this context, the FPZ is represented by a circle centered on a fracture tip [14] (Figure 1e), and r c is the radius of a circular FPZ. A large r c would indicate that micro-cracks are distributed within a large FPZ in front of a preexisting crack tip. It would also suggest a considerable number of microcracks within the FPZ because a preexisting crack will not propagate until the micro-crack density is high enough to reach a critical level [16]. Therefore, r c may characterize the maximum number of micro-cracks generated in the preparation stage for macroscopic crack propagation. In other words, r c can be used to indicate the ability of a rock to tolerate micro-cracks before crack unstable propagation. For this reason, we refer to r c as the crack tolerance. The crack morphology may also depend on the crack tolerance because a large r c would indicate an extensive distribution of micro-cracks, which would likely result in irregular and branch cracks.
The concepts of the FPZ and r c derive from the propagation process of a single crack with specific boundary conditions. Nevertheless, this process represents the inherent mechanical rule of crack generation within rocks because each crack started as an FPZ. Based on this understanding, the crack tolerance may reflect the potential of a given rock stratum to generate abundant cracks. Recent studies [17,18] have shown that rock specimens with a larger FPZ produce more fragments, which suggests greater fracability and provides evidence supporting our hypothesis. The maximum FPZ radius has exhibited dependence on the structure [19,20] and loading conditions [21] of rock. Thus, we conducted experiments and numerical simulations to analyze their effects on crack tolerance and demonstrate its rationality as an evaluation index.

Specimens
We used marble, shale, and sandstone collected from Xishan, Beijing for experiments because marble was observed in some geothermal reservoirs, and shale and sandstone are representative lithologies comprising oil and gas reservoirs. The marble was divided into types A and J (Table 1): marble A totally constituted by calcite had a greater average grain size and was more heterogeneous as defined by Han et al. [22], and marble J mainly consisting of dolomite had an equigranular texture. The microscopy observation and X-ray diffraction (XRD) analysis showed that the shale with fine grains consisted of quartz (55.4%), plagioclase (6.2%), and clay minerals (38.4%, brown grains in Table 1). The quartz and clay minerals were alternatively layered. The sandstone consisted of quartz (69.5%), plagioclase (22.1%), and potassium feldspar (8.4%), and these xenomorphic grains have similar sizes (~2 mm). Most plagioclase grains experienced sericitization. where KIC is the tensile fracture toughness, and rc is the maximum FPZ size. In this context, the FPZ is represented by a circle centered on a fracture tip [14] (Figure 1e), and rc is the radius of a circular FPZ. A large rc would indicate that micro-cracks are distributed within a large FPZ in front of a preexisting crack tip. It would also suggest a considerable number of micro-cracks within the FPZ because a preexisting crack will not propagate until the micro-crack density is high enough to reach a critical level [16]. Therefore, rc may characterize the maximum number of micro-cracks generated in the preparation stage for macroscopic crack propagation. In other words, rc can be used to indicate the ability of a rock to tolerate micro-cracks before crack unstable propagation. For this reason, we refer to rc as the crack tolerance. The crack morphology may also depend on the crack tolerance because a large rc would indicate an extensive distribution of micro-cracks, which would likely result in irregular and branch cracks.
The concepts of the FPZ and rc derive from the propagation process of a single crack with specific boundary conditions. Nevertheless, this process represents the inherent mechanical rule of crack generation within rocks because each crack started as an FPZ. Based on this understanding, the crack tolerance may reflect the potential of a given rock stratum to generate abundant cracks. Recent studies [17,18] have shown that rock specimens with a larger FPZ produce more fragments, which suggests greater fracability and provides evidence supporting our hypothesis. The maximum FPZ radius has exhibited dependence on the structure [19,20] and loading conditions [21] of rock. Thus, we conducted experiments and numerical simulations to analyze their effects on crack tolerance and demonstrate its rationality as an evaluation index.

Specimens
We used marble, shale, and sandstone collected from Xishan, Beijing for experiments because marble was observed in some geothermal reservoirs, and shale and sandstone are representative lithologies comprising oil and gas reservoirs. The marble was divided into types A and J (Table 1): marble A totally constituted by calcite had a greater average grain size and was more heterogeneous as defined by Han et al. [22], and marble J mainly consisting of dolomite had an equigranular texture. The microscopy observation and X-ray diffraction (XRD) analysis showed that the shale with fine grains consisted of quartz (55.4%), plagioclase (6.2%), and clay minerals (38.4%, brown grains in Table 1). The quartz and clay minerals were alternatively layered. The sandstone consisted of quartz (69.5%), plagioclase (22.1%), and potassium feldspar (8.4%), and these xenomorphic grains have similar sizes (~2 mm). Most plagioclase grains experienced sericitization. where KIC is the tensile fracture toughness, and rc is the maximum FPZ size. In this context, the FPZ is represented by a circle centered on a fracture tip [14] (Figure 1e), and rc is the radius of a circular FPZ. A large rc would indicate that micro-cracks are distributed within a large FPZ in front of a preexisting crack tip. It would also suggest a considerable number of micro-cracks within the FPZ because a preexisting crack will not propagate until the micro-crack density is high enough to reach a critical level [16]. Therefore, rc may characterize the maximum number of micro-cracks generated in the preparation stage for macroscopic crack propagation. In other words, rc can be used to indicate the ability of a rock to tolerate micro-cracks before crack unstable propagation. For this reason, we refer to rc as the crack tolerance. The crack morphology may also depend on the crack tolerance because a large rc would indicate an extensive distribution of micro-cracks, which would likely result in irregular and branch cracks.
The concepts of the FPZ and rc derive from the propagation process of a single crack with specific boundary conditions. Nevertheless, this process represents the inherent mechanical rule of crack generation within rocks because each crack started as an FPZ. Based on this understanding, the crack tolerance may reflect the potential of a given rock stratum to generate abundant cracks. Recent studies [17,18] have shown that rock specimens with a larger FPZ produce more fragments, which suggests greater fracability and provides evidence supporting our hypothesis. The maximum FPZ radius has exhibited dependence on the structure [19,20] and loading conditions [21] of rock. Thus, we conducted experiments and numerical simulations to analyze their effects on crack tolerance and demonstrate its rationality as an evaluation index.

Specimens
We used marble, shale, and sandstone collected from Xishan, Beijing for experiments because marble was observed in some geothermal reservoirs, and shale and sandstone are representative lithologies comprising oil and gas reservoirs. The marble was divided into types A and J (Table 1): marble A totally constituted by calcite had a greater average grain size and was more heterogeneous as defined by Han et al. [22], and marble J mainly consisting of dolomite had an equigranular texture. The microscopy observation and X-ray diffraction (XRD) analysis showed that the shale with fine grains consisted of quartz (55.4%), plagioclase (6.2%), and clay minerals (38.4%, brown grains in Table 1). The quartz and clay minerals were alternatively layered. The sandstone consisted of quartz (69.5%), plagioclase (22.1%), and potassium feldspar (8.4%), and these xenomorphic grains have similar sizes (~2 mm). Most plagioclase grains experienced sericitization.

Experimental Methodology
The cracked chevron notched Brazilian disk (CCNBD) test involves the formation of FPZs at the two tips of a prefabricated notched crack, which is analogous to a natural crack. Therefore, the CCNBD test is applicable to evaluating crack tolerance. According to Equation (3), quantifying the crack tolerance requires determining the tensile fracture toughness KIC and tensile strength σt, which are measured by the CCNBD and Brazilian disk (BD) tests, respectively, as recommended by the International Society for Rock Mechanics (ISRM) and American Society of Testing Materials (ASTM).
The notched crack of each CCNBD specimen was created by a 1 mm thick circular diamond saw. To ensure cutting accuracy, the expected locations of the circular center and the initial and final chevron notched cracks were marked on each disk. We measured the actual values of the parameters shown in Figure 2a,b and confirmed that the dimensionless parameters α1 and αB of all CCNBD specimens were within the valid range ( Figure  2c). The method reported by Fowell et al. [23] was used to calculate the KIC: where Pmax is the peak applied axial load in the CCNBD test and Y * min is the critical dimensionless stress intensity value. This is determined by where u and v are geometric constants that are determined by α0 and αB as reported by Fowell et al. [23].

Experimental Methodology
The cracked chevron notched Brazilian disk (CCNBD) test involves the formation of FPZs at the two tips of a prefabricated notched crack, which is analogous to a natural crack. Therefore, the CCNBD test is applicable to evaluating crack tolerance. According to Equation (3), quantifying the crack tolerance requires determining the tensile fracture toughness KIC and tensile strength σt, which are measured by the CCNBD and Brazilian disk (BD) tests, respectively, as recommended by the International Society for Rock Mechanics (ISRM) and American Society of Testing Materials (ASTM).
The notched crack of each CCNBD specimen was created by a 1 mm thick circular diamond saw. To ensure cutting accuracy, the expected locations of the circular center and the initial and final chevron notched cracks were marked on each disk. We measured the actual values of the parameters shown in Figure 2a,b and confirmed that the dimensionless parameters α1 and αB of all CCNBD specimens were within the valid range ( Figure  2c). The method reported by Fowell et al. [23] was used to calculate the KIC: where Pmax is the peak applied axial load in the CCNBD test and Y * min is the critical dimensionless stress intensity value. This is determined by where u and v are geometric constants that are determined by α0 and αB as reported by Fowell et al. [23].

Experimental Methodology
The cracked chevron notched Brazilian disk (CCNBD) test involves the formation of FPZs at the two tips of a prefabricated notched crack, which is analogous to a natural crack. Therefore, the CCNBD test is applicable to evaluating crack tolerance. According to Equation (3), quantifying the crack tolerance requires determining the tensile fracture toughness K IC and tensile strength σ t , which are measured by the CCNBD and Brazilian disk (BD) tests, respectively, as recommended by the International Society for Rock Mechanics (ISRM) and American Society of Testing Materials (ASTM).
The notched crack of each CCNBD specimen was created by a 1 mm thick circular diamond saw. To ensure cutting accuracy, the expected locations of the circular center and the initial and final chevron notched cracks were marked on each disk. We measured the actual values of the parameters shown in Figure 2a,b and confirmed that the dimensionless parameters α 1 and α B of all CCNBD specimens were within the valid range ( Figure 2c). The method reported by Fowell et al. [23] was used to calculate the K IC : where P max is the peak applied axial load in the CCNBD test and Y * min is the critical dimensionless stress intensity value. This is determined by where u and v are geometric constants that are determined by α 0 and α B as reported by Fowell et al. [23]. The thickness (B') and diameter (D') of the BD specimens were set identical to those of the CCNBD specimens to eliminate the size effect on the calculated crack tolerance. The B'-to-D' ratio was within the range recommended by the ASTM of 0.2-0.75 [24]. The σ t was calculated as follows: where P' max is the peak applied axial load in the BD test. The thickness (B') and diameter (D') of the BD specimens were set identical to those of the CCNBD specimens to eliminate the size effect on the calculated crack tolerance. The B'-to-D' ratio was within the range recommended by the ASTM of 0.2-0.75 [24]. The σt was calculated as follows: where P'max is the peak applied axial load in the BD test.
Each CCNBD or BD test (Figure 3a,b) was performed at a constant displacement rate of 0.06 mm/min by an MTS servo-control testing machine (series CMT) with a maximum loading force of 100 kN. This machine is equipped with an SNAS GDS-300 environmental chamber controlled by a WK650 controller (Figure 3c,d). These apparatuses permit environmental temperatures within the chamber up to 200 °C by electrical heaters (Figure 3b). To investigate the effect of temperature, several sandstone specimens were placed in the chamber at 75 or 125 °C for 1 h before the tests began, so that the notched crack propagated within rocks under higher temperatures. Other tests were performed at room temperature (~25 °C). The bedding planes of the shale specimens were set perpendicular (horizontal) or parallel (vertical) to the notched cracks to analyze the effect of the bedding orientation. Each CCNBD or BD test (Figure 3a,b) was performed at a constant displacement rate of 0.06 mm/min by an MTS servo-control testing machine (series CMT) with a maximum loading force of 100 kN. This machine is equipped with an SNAS GDS-300 environmental chamber controlled by a WK650 controller (Figure 3c,d). These apparatuses permit environmental temperatures within the chamber up to 200 • C by electrical heaters (Figure 3b). To investigate the effect of temperature, several sandstone specimens were placed in the chamber at 75 or 125 • C for 1 h before the tests began, so that the notched crack propagated within rocks under higher temperatures. Other tests were performed at room temperature (~25 • C). The bedding planes of the shale specimens were set perpendicular (horizontal) or parallel (vertical) to the notched cracks to analyze the effect of the bedding orientation. BD tests were conducted on at least three specimens in parallel with the same lithology, bedding orientation, and temperature, and the average strength was taken as the tensile strength for the corresponding set of conditions. The KIC of each CCNBD specimen and the above average σt were used in Equation (3) to calculate the crack tolerance. BD tests were conducted on at least three specimens in parallel with the same lithology, bedding orientation, and temperature, and the average strength was taken as the tensile strength for the corresponding set of conditions. The K IC of each CCNBD specimen and the above average σ t were used in Equation (3) to calculate the crack tolerance.

Experimental Results
In the CCNBD test, the marble A specimens with stronger heterogeneity had greater tensile strength, fracture toughness, and crack tolerance than the marble J specimens, with relatively homogeneous small grains (Figure 4a). White patches indicating FPZs [25] appeared in front of the notched crack tips (Figure 5a) as the peak loads of the marble A specimens were approached. The patches corresponded to the sparkling areas on the rupture surface (Figure 5b), which may imply breaking cleavages of grains. However, such patches were not observed for the marble J specimens (Figure 5c,d), and neither were the discernible sparkling areas. Furthermore, the main cracks in the marble A specimens had branches causing more fragments (Figure 5b) while the crack in the marble J specimens propagated along a straight path (Figure 5c). These phenomena suggest that a more heterogeneous grain size corresponds to a larger crack tolerance and thus a stronger ability for crack generation.     The mean crack tolerance of the shale specimens was less with a vertical bedding orientation than with a horizontal orientation (Figure 4b). The tensile strength and fracture toughness displayed similar variation trends with bedding orientation. Similar results can be acquired based on the data from Wang [26]. With a vertical orientation, the main crack of the specimen propagated along the bedding planes (Figure 6a), which generated a smooth rupture surface (Figure 6b). In contrast, with a horizontal orientation, the main crack spanned across bedding planes, and the path with steps was more irregular (Figure 6c,d). This is because the main crack was offset or even bifurcated when it encountered a bedding plane. The branch cracks were captured by bedding planes and then propagated along them, thereby their morphologies were smooth. The mean crack tolerance of the shale specimens was less with a vertical bedding orientation than with a horizontal orientation (Figure 4b). The tensile strength and fracture toughness displayed similar variation trends with bedding orientation. Similar results can be acquired based on the data from Wang [26]. With a vertical orientation, the main crack of the specimen propagated along the bedding planes (Figure 6a), which generated a smooth rupture surface (Figure 6b). In contrast, with a horizontal orientation, the main crack spanned across bedding planes, and the path with steps was more irregular (Figure 6c,d). This is because the main crack was offset or even bifurcated when it encountered a bedding plane. The branch cracks were captured by bedding planes and then propagated along them, thereby their morphologies were smooth. The crack tolerance value of the sandstone specimens consistently declined as the environmental temperature rises from 25 °C to 125 °C, while the tensile strength and fracture toughness exhibited V-shaped trends within this temperature range (Figure 4c). It is difficult to identify changes in crack morphology with the rising temperature with the naked eye (Figure 7a,c,e). According to the edge of their rupture surface, we speculated The crack tolerance value of the sandstone specimens consistently declined as the environmental temperature rises from 25 • C to 125 • C, while the tensile strength and fracture toughness exhibited V-shaped trends within this temperature range (Figure 4c). It is difficult to identify changes in crack morphology with the rising temperature with the naked eye (Figure 7a,c,e). According to the edge of their rupture surface, we speculated that the main cracks in the specimens at 125 • C may propagate along less curved paths than the specimens at lower temperatures did (Figure 7b,d,f). The variations of the crack tolerance value and crack morphologies imply that high temperatures possibly reduce rock fracability.

Particle Flow Code
To test the rationality of the crack tolerance as an evaluation index for rock fracability, we adopt particle flow code in two dimensions (PFC 2D ), which is widely used for discrete element method (DEM). Rock was modeled as a dense packing of non-uniform-sized and inter-bonded circular particles using this method, and its mechanical behavior relied on the microscale properties and constitutive relations of the bonded contacts between the particles. Following Newton's laws of motion, the force acting at each contact were updated with the particle movements during the simulation process, and the breakage of bonds representing crack generation [28] occurred when a component of the contact force satisfied a certain criterion.
We used the experimental results for marble A and J as examples for the DEM simulation because the marble contained polygonal minerals >1 mm in size, which allowed us to implement a grain-based model (GBM) with a polygon-tessellation grain boundaries [29]. Such a model takes the mineral grain texture into account (Figure 8a), making a simulation more vivid. The modeling method for GBM of marble refers to [27]. Soft-bonded [30] and smooth-joint models [31] were employed to express bonded and unbonded behaviors characteristic of intra-grain and inter-grain contacts (Figure 8b-d), respectively. The crack tolerance preliminarily showed an ability to address the aforementioned problems of previous fracability indices. Firstly, this index has a firm physical basis derived from the FPZ size, representing the nonlinearity of deformation due to micro-crack generation before the macroscopic propagation of the crack. Additionally, the formation of the FPZ is the inherent mechanical behavior of heterogeneous rock materials, and the FPZ size highly depends on the degree of heterogeneity as the previous [27] and present experiment results revealed. From the aspect of practicability, the crack tolerance value can be determined conveniently in the laboratory because BD and CCNBD tests are very common rock mechanical tests, and the small-size specimens they use can be easily obtained from cores. Finally, this index may characterize the effects of structure and loading conditions on fracability to an extent, as the tests on the shale and sandstone showed.

Particle Flow Code
To test the rationality of the crack tolerance as an evaluation index for rock fracability, we adopt particle flow code in two dimensions (PFC 2D ), which is widely used for discrete element method (DEM). Rock was modeled as a dense packing of non-uniform-sized and inter-bonded circular particles using this method, and its mechanical behavior relied on the microscale properties and constitutive relations of the bonded contacts between the particles. Following Newton's laws of motion, the force acting at each contact were updated with the particle movements during the simulation process, and the breakage of bonds representing crack generation [28] occurred when a component of the contact force satisfied a certain criterion.
We used the experimental results for marble A and J as examples for the DEM simulation because the marble contained polygonal minerals >1 mm in size, which allowed us to implement a grain-based model (GBM) with a polygon-tessellation grain boundaries [29]. Such a model takes the mineral grain texture into account (Figure 8a), making a simulation more vivid. The modeling method for GBM of marble refers to [27]. Soft-bonded [30] and smooth-joint models [31] were employed to express bonded and unbonded behaviors characteristic of intra-grain and inter-grain contacts (Figure 8b-d), respectively.

Model Setup and Parameter Calibration
Based on the grain size distributions in Table 1, four circular 75 mm GBMs were created representing the marble A and J specimens in the BD and CCNBD tests. Each model comprised ~20,000 circular basic particles with a 0.2-0.3 mm radius. Since even a single mineral crystal is anisotropic along different atomic lattices [32], we set the strength and deformation parameters of the soft-bonded contacts to follow the Weibull distribution, and the shape parameter representing heterogeneity was set to 3 and 5 for marble A and J, respectively. A small value for the shape parameter indicates strong heterogeneity.
Before conducting the simulation of the CCNBD tests, the microscale parameters of the particles and contacts required iteratively calibrating through trial and error referring to the BD test results and the previous work [27,33,34]. The GBMs of the BD specimens were positioned between two stiff walls representing the loading end and platform of a compression machine, and the walls moved toward each other at the same constant velocity to result in a quasi-static loading rate. The calibration completes until the simulated load-displacement curves and crack morphology fit well with the observations in the BD tests ( Figure 9). The calibrated microscale parameters ( Table 2) were used to simulate the marble specimens in the CCNBD test.

Model Setup and Parameter Calibration
Based on the grain size distributions in Table 1, four circular 75 mm GBMs were created representing the marble A and J specimens in the BD and CCNBD tests. Each model comprised~20,000 circular basic particles with a 0.2-0.3 mm radius. Since even a single mineral crystal is anisotropic along different atomic lattices [32], we set the strength and deformation parameters of the soft-bonded contacts to follow the Weibull distribution, and the shape parameter representing heterogeneity was set to 3 and 5 for marble A and J, respectively. A small value for the shape parameter indicates strong heterogeneity.
Before conducting the simulation of the CCNBD tests, the microscale parameters of the particles and contacts required iteratively calibrating through trial and error referring to the BD test results and the previous work [27,33,34]. The GBMs of the BD specimens were positioned between two stiff walls representing the loading end and platform of a compression machine, and the walls moved toward each other at the same constant velocity to result in a quasi-static loading rate. The calibration completes until the simulated loaddisplacement curves and crack morphology fit well with the observations in the BD tests ( Figure 9). The calibrated microscale parameters (

Simulation Results
In the numerical simulations of the CCNBD test, when the applied load reached a certain level, micro-cracking was initiated near the notched crack tips of the specimens (Figures 10a and 11a). At the peak loads (Pmax) of the marble A and J specimens, the microcracks around the crack tips tended to coalesce to form new macroscopic cracks (Figures  10b and 11b). After that, the notched crack propagated dramatically, which caused a rapid post-peak drop in the applied load and the specimen to rupture (Figures 10c and 11c). Therefore, the preparation stage for dramatic propagation of a notched crack can be defined as from the initiation of micro-cracking to the reaching of the peak load, during which micro-cracks generate to develop the FPZ. The FPZ is the area near the crack tip with a dense micro-crack distribution when the peak load is reached that stays in the critical state of macroscopic rupture. As mentioned in Section 2, the crack tolerance is characterized by the size of the FPZ.

Simulation Results
In the numerical simulations of the CCNBD test, when the applied load reached a certain level, micro-cracking was initiated near the notched crack tips of the specimens (Figures 10a and 11a). At the peak loads (P max ) of the marble A and J specimens, the micro-cracks around the crack tips tended to coalesce to form new macroscopic cracks (Figures 10b and 11b). After that, the notched crack propagated dramatically, which caused a rapid post-peak drop in the applied load and the specimen to rupture (Figures 10c and 11c). Therefore, the preparation stage for dramatic propagation of a notched crack can be defined as from the initiation of micro-cracking to the reaching of the peak load, during which micro-cracks generate to develop the FPZ. The FPZ is the area near the crack tip with a dense micro-crack distribution when the peak load is reached that stays in the critical state of macroscopic rupture. As mentioned in Section 2, the crack tolerance is characterized by the size of the FPZ.  As the simulation results (Figures 10b and 11b) showed, the micro-crack density at the vicinity of notched crack tip was especially high in the whole specimen, owing to the nonlinear deformation brought by stress concentration. The FPZ was assumed as a tipcentered circle with radius of rc [14] that covered the area with high micro-crack density. Therefore, with increasing distance (radius) Rf from the tip and the diminishing intensity of the stress concentration, the deformation transitions from nonlinear inside the FPZ to quasilinear outside the FPZ, and thus the micro-crack density outside the FPZ declined to the background density of the rock [35].
To simplify the analysis, we assumed that the distributions of micro-crack inside and outside the FPZ are uniform but have different density. On the basis of this, the total micro-crack number N within a certain circular statistical range with Rf radius can be formulated as where d and d0 are the average micro-crack density inside and outside the FPZ, respectively. N displayed a positively correlation with Rf; however, the curves of N deflected when Rf increased to rc that defined the boundary of the FPZ (Figure 12). This is because the micro-crack density d inside the FPZ can be up to ~15 times as great as the background density d0 [36]; the increasing rate of N will decelerate once the statistical range extend outside the FPZ. Such a deflection became more identifiable with the increasing ratio of d/d0. Thus, the crack tolerance of the specimens can be determined by the radius corresponding to the deflection point.  As the simulation results (Figures 10b and 11b) showed, the micro-crack density at the vicinity of notched crack tip was especially high in the whole specimen, owing to the nonlinear deformation brought by stress concentration. The FPZ was assumed as a tipcentered circle with radius of rc [14] that covered the area with high micro-crack density. Therefore, with increasing distance (radius) Rf from the tip and the diminishing intensity of the stress concentration, the deformation transitions from nonlinear inside the FPZ to quasilinear outside the FPZ, and thus the micro-crack density outside the FPZ declined to the background density of the rock [35].
To simplify the analysis, we assumed that the distributions of micro-crack inside and outside the FPZ are uniform but have different density. On the basis of this, the total micro-crack number N within a certain circular statistical range with Rf radius can be formulated as where d and d0 are the average micro-crack density inside and outside the FPZ, respectively. N displayed a positively correlation with Rf; however, the curves of N deflected when Rf increased to rc that defined the boundary of the FPZ (Figure 12). This is because the micro-crack density d inside the FPZ can be up to ~15 times as great as the background density d0 [36]; the increasing rate of N will decelerate once the statistical range extend outside the FPZ. Such a deflection became more identifiable with the increasing ratio of d/d0. Thus, the crack tolerance of the specimens can be determined by the radius corresponding to the deflection point. As the simulation results (Figures 10b and 11b) showed, the micro-crack density at the vicinity of notched crack tip was especially high in the whole specimen, owing to the nonlinear deformation brought by stress concentration. The FPZ was assumed as a tip-centered circle with radius of r c [14] that covered the area with high micro-crack density. Therefore, with increasing distance (radius) R f from the tip and the diminishing intensity of the stress concentration, the deformation transitions from nonlinear inside the FPZ to quasilinear outside the FPZ, and thus the micro-crack density outside the FPZ declined to the background density of the rock [35].
To simplify the analysis, we assumed that the distributions of micro-crack inside and outside the FPZ are uniform but have different density. On the basis of this, the total micro-crack number N within a certain circular statistical range with R f radius can be formulated as where d and d 0 are the average micro-crack density inside and outside the FPZ, respectively. N displayed a positively correlation with R f ; however, the curves of N deflected when R f increased to r c that defined the boundary of the FPZ (Figure 12). This is because the microcrack density d inside the FPZ can be up to~15 times as great as the background density d 0 [36]; the increasing rate of N will decelerate once the statistical range extend outside the FPZ. Such a deflection became more identifiable with the increasing ratio of d/d 0 . Thus, the crack tolerance of the specimens can be determined by the radius corresponding to the deflection point. As the above calculation predicts, the micro-crack number N in Figures 10b and 11b for the marble A and J specimens increased with Rf, and the N-Rf curves deflected at radii of 7 and 4 mm, respectively (Figure 13), which were closed to the mean values of the crack tolerance of A (~8 mm) and J (~3 mm) measured by experiments. These results demonstrate that calculating the crack tolerance using the tensile fracture toughness and average tensile strength in Equation (3) leads to reliable results. The simulations also showed that marble A had a greater crack tolerance than marble J, and the FPZ of the former contained more micro-cracks than that of the latter preceding specimen rupture (Figure 13c,d). Correspondingly, the GBM of marble A generates 1473 micro-cracks in the whole loading process, more than that of marble J (1404). These results showed that a greater crack tolerance can represent a stronger ability to generate micro-cracks.  As the above calculation predicts, the micro-crack number N in Figures 10b and 11b for the marble A and J specimens increased with R f , and the N-R f curves deflected at radii of 7 and 4 mm, respectively ( Figure 13), which were closed to the mean values of the crack tolerance of A (~8 mm) and J (~3 mm) measured by experiments. These results demonstrate that calculating the crack tolerance using the tensile fracture toughness and average tensile strength in Equation (3) leads to reliable results. The simulations also showed that marble A had a greater crack tolerance than marble J, and the FPZ of the former contained more micro-cracks than that of the latter preceding specimen rupture (Figure 13c,d). Correspondingly, the GBM of marble A generates 1473 micro-cracks in the whole loading process, more than that of marble J (1404). These results showed that a greater crack tolerance can represent a stronger ability to generate micro-cracks. As the above calculation predicts, the micro-crack number N in Figures 10b and 11b for the marble A and J specimens increased with Rf, and the N-Rf curves deflected at radii of 7 and 4 mm, respectively ( Figure 13), which were closed to the mean values of the crack tolerance of A (~8 mm) and J (~3 mm) measured by experiments. These results demonstrate that calculating the crack tolerance using the tensile fracture toughness and average tensile strength in Equation (3) leads to reliable results. The simulations also showed that marble A had a greater crack tolerance than marble J, and the FPZ of the former contained more micro-cracks than that of the latter preceding specimen rupture (Figure 13c,d). Correspondingly, the GBM of marble A generates 1473 micro-cracks in the whole loading process, more than that of marble J (1404). These results showed that a greater crack tolerance can represent a stronger ability to generate micro-cracks.  The strike angle of micro-cracks ( Figure 14) for the marble J mainly distributed in the range of 70-100 • , which was narrower than that, 60-110 • , for the marble A. This result suggests that more micro-cracks deviated from the loading direction (90 • ) in the marble A. The coalescence of such micro-cracks with various strike angles resulted in macroscopic cracks that propagated along irregular even branched paths and radiated from the notched crack tips (Figures 5 and 10). Otherwise, the macro-cracks will develop primarily parallel to the loading direction, and thus their morphologies were less complex, as the marble J specimen showed. Therefore, these observations confirmed the assumption that fracability can be characterized by the crack tolerance.
Materials 2022, 15, x FOR PEER REVIEW 15 The strike angle of micro-cracks ( Figure 14) for the marble J mainly distributed i range of 70°-100°, which was narrower than that, 60°-110°, for the marble A. This r suggests that more micro-cracks deviated from the loading direction (90°) in the m A. The coalescence of such micro-cracks with various strike angles resulted in m scopic cracks that propagated along irregular even branched paths and radiated from notched crack tips (Figures 5 and 10). Otherwise, the macro-cracks will develop prim parallel to the loading direction, and thus their morphologies were less complex, a marble J specimen showed. Therefore, these observations confirmed the assumption fracability can be characterized by the crack tolerance.

Effect of Grain Size on Crack Tolerance
The grain size greatly influences the cracking behavior of rocks [37]. Regarding consisting of grains with various size, the grain size greatly differs across parts of the specimen, so does the microscopic strength, which enhances the rock heterogeneity. is why the heterogeneity index (shape parameter for the Weibull distribution) of m A was set as smaller than that for marble J in the GBMs.
The microscopic strength in different parts of a strongly heterogeneous rock s men can distribute in a wide range, so a small increment in the stress near a crack tip easily cause micro-cracking within such rocks. Therefore, the initiation of micro-crac was earlier, i.e., corresponding to a smaller ratio between the applied load and the load, in marble A than in marble J (Figures 10 and 11). However, cracking can als arrested easily because it probably encounters stronger local parts soon. Therefore rupture of strongly heterogeneous specimens will not occur until there are adequat cro-cracks to create FPZs and macroscopic cracks. In summary, strong heteroge strengthens the crack tolerance of rocks.
For rock specimens with a homogeneous grain size distribution, its microsc strength in different parts can be generally closed to a certain level. Thus, only a few cro-cracks arise before the stress near the notched crack tips reaches that strength l Once the strength is reached, the cracks propagate dramatically, which causes a r rupture. Macroscopic cracks spanning across the specimen form nearly instantaneo following the generation of a small FPZ. For these reasons, relatively homogeneous ble J had a smaller crack tolerance than marble A.

Effects of the Bedding Orientation and Environmental Temperature
With a vertical bedding orientation, micro-cracking naturally initiates within ding planes in front of crack tips and propagates along them because the tensile stre

Effect of Grain Size on Crack Tolerance
The grain size greatly influences the cracking behavior of rocks [37]. Regarding rock consisting of grains with various size, the grain size greatly differs across parts of the rock specimen, so does the microscopic strength, which enhances the rock heterogeneity. This is why the heterogeneity index (shape parameter for the Weibull distribution) of marble A was set as smaller than that for marble J in the GBMs.
The microscopic strength in different parts of a strongly heterogeneous rock specimen can distribute in a wide range, so a small increment in the stress near a crack tip can easily cause micro-cracking within such rocks. Therefore, the initiation of micro-cracking was earlier, i.e., corresponding to a smaller ratio between the applied load and the peak load, in marble A than in marble J (Figures 10 and 11). However, cracking can also be arrested easily because it probably encounters stronger local parts soon. Therefore, the rupture of strongly heterogeneous specimens will not occur until there are adequate micro-cracks to create FPZs and macroscopic cracks. In summary, strong heterogeneity strengthens the crack tolerance of rocks.
For rock specimens with a homogeneous grain size distribution, its microscopic strength in different parts can be generally closed to a certain level. Thus, only a few micro-cracks arise before the stress near the notched crack tips reaches that strength level. Once the strength is reached, the cracks propagate dramatically, which causes a rapid rupture. Macroscopic cracks spanning across the specimen form nearly instantaneously following the generation of a small FPZ. For these reasons, relatively homogeneous marble J had a smaller crack tolerance than marble A.

Effects of the Bedding Orientation and Environmental Temperature
With a vertical bedding orientation, micro-cracking naturally initiates within bedding planes in front of crack tips and propagates along them because the tensile strength of shale bedding planes is usually much smaller than that of layers between planes. Therefore, the tensile strength and fracture toughness were relatively low. Since the micro-cracks are limited to thin bedding planes, the corresponding crack tolerance is also small. With a horizontal bedding orientation, layers comprising various mineral grains make the distribution of the microscopic tensile strength near each tip more heterogeneous, so microcracking may be scattered among the layers. The behavior of micro-cracks can be complex when they cross bedding planes: cracks may branch along bedding planes and result in curved macroscopic cracks. Although hydraulic fracturing involves with many factors affecting the interaction between cracks and bedding planes [38], it is generally recognized that a crack propagating along a bedding plane is the most unfavorable situation for generating complex fracture networks [39,40].
Thermal treatment of sandstone leads to dehydration and the thermal expansion of minerals, which promotes the brittle-ductile transition of minerals [41]. Dehydration occurs at~100 • C, which is when absorbed water escapes from the mineral surface [42], and enhances the friction and bonding strength between minerals. Thermal expansion takes effect when the temperature exceeds 100 • C and closes preexisting micro-cracks in rocks [43], which enhances the tensile strength of mineral grains and boundaries within a certain temperature range. Owing to these effects arising exceeding 100 • C, the sandstone at 125 • C had higher tensile strength and fracture toughness than that at 75 • C (Figure 4c). However, the dehydration and thermal expansion of minerals reduce the structural heterogeneity and rock fracability. In addition, fracability also weakens when sandstone transitions from brittle to ductile [44]. These observations suggest that increasing the temperature from 25 • C to 125 • C should be unfavorable for crack generation in sandstone. Such a negative effect of temperature was also observed in Longmaxi shale [45], a commercial shale gas reservoir in Chongqing, China. Considering that the downhole temperature, especially in deep and geothermal wells, can be much higher than the surface temperature, more attention should be given to the effect of temperature over a broader range on rock fracability.

Implications in Hydraulic Fracturing
The FPZ indicates nonlinear deformation (i.e., micro-crack generation) within a rock, which originates from rock heterogeneity. Nonlinear deformation diminishes and transitions into linear deformation with decreasing rock heterogeneity, which would reduce the FPZ size (i.e., crack tolerance). Complex cracks barely form in the absence of micro-cracks and their interactions. Because the formation of the FPZ is intrinsic to heterogeneous rock, the association between the FPZ size and heterogeneity is the physical basis for the crack tolerance as an evaluation index of fracability. The crack tolerance can reflect the effects of structural and environmental factors because they influence rock heterogeneity [46,47].
A large crack tolerance indicates many micro-cracks within the FPZ, which would tend to cause an irregular morphology of macroscopic cracks and wider zones of micro-cracks along both sides of cracks. These characteristics allow more cracks and pores in the rocks to connect with the main cracks during crack propagation driven by fluid. This increases the volume of cracks, which enhances the fluid conductivity of rock to exploit oil, gas, and geothermal resources.
The present study mainly provided laboratory observations on the effects of three factors to support the rationality of crack tolerance in reservoir fracability evaluation. However, the complexity of crack networks is dependent on various factors. Further investigations that consider other environmental and structural effects (e.g., magnitude and direction of crustal stress) are required to test this fracability index. On the other hand, the reliability and practical applicability of this index need to be further tested using rock specimens collected from reservoirs being exploited.

Conclusions
Aiming to the challenge of lacking a reliable index of fracability evaluation for hydraulic fracturing, we suggest that the crack tolerance, i.e., the maximum radius of the FPZ, may be used to characterize the fracability.
Crack tolerance originates from the unique heterogeneity of rock and inherent rules of crack generation, and thus it has a clear physical meaning and firm mechanical basis.
This index can be conveniently quantified in the laboratory using BD and CCNBD tests. We showed that the crack tolerance is positively correlated with the grain size heterogeneity and is negatively correlated with the environmental temperature (25-125 • C). The crack tolerance of shale specimens was greater with a horizontal bedding orientation than with a vertical orientation. Thus, crack tolerance can well characterize the effects of certain rock properties and environmental factors.
In summary, crack tolerance is promising in serving as a reliable evaluation index of rock fracability in terms of rationality and practicability, which has significant engineering implications for efficient hydraulic fracturing.  Data Availability Statement: All the data required to evaluate the conclusions of this study are present in the paper. The authors will provide additional data related to this paper upon request.

Conflicts of Interest:
The authors declare no conflict of interest.