Mutual interference of layer plane and natural fracture in the failure behavior of shale and the mechanism investigation

Shale contains a certain amount of natural fractures, which affects the mechanical properties of shale. In this paper, a bonded-particle model in particle flow code (PFC) is established to simulate the failure process of layered shale under Brazilian tests, under the complex relationship between layer plane and natural fracture. First, a shale model without natural fractures is verified against the experimental results. Then, a natural fracture is embedded in the shale model, where the outcomes indicate that the layer plane angle (marked as α) and the angle (marked as β) of embedded fracture prominently interfere the failure strength anisotropy and fracture pattern. Finally, sensitivity evaluations suggest that variable tensile/cohesion strength has a changeable influence on failure mechanism of shale, even for same α or/and β. To serve this work, the stimulated fractures are categorized into two patterns based on whether they relate to natural fracture or not. Meanwhile, four damage modes and the number of microcracks during the loading process are recognized quantitatively to study the mechanism of shale failure behavior. Considering the failure mechanism determines the outcome of hydraulic fracturing in shale, this work is supposed to provide a significant implication in theory for the engineering operation.


Introduction
Hydraulic fracturing is an effective technique to stimulate shale reservoirs, which enables the boom of shale gas/oil development over the world Yin et al. 2017;Liu et al. 2019aLiu et al. , b, 2020Zhou et al. 2019;Zhao et al. 2020a). Shale is usually characterized as layered and contains complex natural fractures (Liu et al. 2016(Liu et al. , 2017. Previous achievements noted that these natural fractures make the mechanical properties of shale more complex, resulting in that the failure patterns of shale mainly depend on the mechanical characteristics of the layer planes and natural fractures as well as the relationship between them (e.g., März et al. 2008;Vervoort et al. 2014;Tan et al. 2015;Feng et al. 2020). Therefore, there is a strong need to better understand the mechanical characteristics and failure process of shale regarding fracturing engineering applications (Li et al. , 2020Zhao et al. 2020b).
Natural fractures are variable to shale with a different geological background (e.g., sedimentary or structural environment), resulting in the diverse properties of natural fractures, such as tension strength, cohesion and friction. Therefore, it is challenging to universally describe the failure mechanisms of shale containing natural fracture in theory. To address this issue, rock mechanics testing methods such as Brazilian tests, compression tests and notched three-pointbending tests were conducted in laboratory to investigate the failure behavior of rocks with flaws, in which these flaws simulate natural fractures (Roy et al. 2017;Na et al. 2017;Yang et al. 2017;Feng et al. 2019). For shale rock, Yang et al. (2015), Zhang et al. (2018a, b) and Nezhad et al. (2018) experimentally studied the effect of the layer-loading angle on the failure mode of Brazilian tests. Heng et al. (2015) conducted direct shear tests to study the mechanical properties of layer planes and the shear strength of shale with different layer orientations. More examples like Arora Edited by Xiu-Qiu Peng and Mishra (2015), Geng et al. (2016) and Al-Maamori et al. (2019) made triaxial compression tests on shale to analyze the variation law of transversely isotropic elastic parameters, aiming to obtain the relationship between the elastic parameters and layer angle. Shi et al. (2019) conducted notched three-point bending tests to elucidate the fracture behavior of anisotropic shale. By experimental studies, these above researches accumulated a lot of knowledge about mechanical characteristic of shale. However, in laboratory, it is not easy to prepare shale samples containing flaws (preexisting fracture) with variable characteristics (e.g., width, length, intersection angle with layer plane, etc.). This shortage of experimental attempts limits researchers to systematically interpret the failure property of layered shale with natural fracture, resulting in that the effect of natural fracture on the shale failure behavior remains unclear and needs to be examined.
As an alternative research method, numerical simulation can be regarded as a supplement to the experiments, since it is able to easily control the variables and gather more operating conditions in a model. According to different means of describing weak planes, existing numerical methods are divided into continuous medium methods (CMMs) and discrete element methods (DEMs). In contrast to the CMM (Cai 2013;Zhu et al. 2019), DEM procedure is able to intuitively embody the weak planes, where the initiation and propagation of fractures can be obtained explicitly without applying complex constitutive laws (Bennett et al. 2015;Park and Min 2015;Zhang et al. 2015;Wang et al. 2017Wang et al. , 2018Zhang et al. 2018a, b;Peng et al. 2018). In addition, considering the small scale (e.g., diameter of 50 mm) of the specimens, the particle DEM, e.g., particle flow code (PFC), is widely used to describe the gradual failure of rocks and the microscopic mechanisms underlying rock deformation behavior (Bahaaddini et al. 2016;Jia et al. 2017). Hence, the PFC is widely introduced to the numerical simulation about shale issues. For example, He and Afolagboye (2018) investigated the influence of the shale layer plane orientation and layer cohesion on the tensile strength and fracture modes under Brazilian tests. Chong et al. (2017) proposed an anisotropic mineral brittleness-based model to study shale properties and found that the ratio of cohesion to tensile strength of smooth joints mainly affects the number of cracks formed. Besides, Yang and Huang (2014), Luo et al. (2018) and Dou et al. (2019) also numerically studied the mechanic behavior using PFC. Unfortunately, in spite that numerical investigations were conducted extensively, of which sparing PFC-based ones have gone to the failure performance of shale containing natural fracture.
In this study, a DEM-based numerical model is established for Brazilian tests. Based on PFC2D, the established layered shale model without fractures (referred to as intact shale) is first validated against laboratory results. Then, a single flaw is embedded into the intact model to simulate the failure behavior of shale with natural fracture. Finally, this work further discusses the sensitivity of layer plane strength parameters (namely tensile/cohesion strength) on the failure strength and fracture pattern of shale containing natural fracture. During the investigation, this work puts forward a novel and concise classification scheme on the fractures stimulated by Brazilian tests. Moreover, this work quantitatively identifies the damage modes (matrix shear, matrix tension, layer shear and layer tension) and the evolution of microcracks during the loading process, aiming to expose the mechanism of shale failure behavior. Basically, this study focuses on the interference of different discontinuities (natural fracture and layer plane) on the failure behavior of shale, which is expected to be helpful in guiding field engineering operation of hydraulic fracture theoretically.

Numerical model establishment and validation
An intact shale model is built for the numerical model validation by compared with the experimental investigation using the Brazilian disk test. This model has a disk shape with a diameter of 50 mm, which is the same as the experimental disks. According to the sensitivity analysis of numerical parameters (Ding et al. 2014), reasonable results can be obtained when the model size is 50 times larger than the average particle radius and the particle radius ratio is close to 1.66. Moreover, shale comprises a series of thin layer planes, and the spacing between the layers generally varies from 0.1 mm to 1 mm (Zhang et al. 2017). To ensure the efficiency of the calculation and the accuracy of the results, the distance between the layers of the numerical model is set to 1 mm. The total number of particles is 14,661, with a uniform distribution of particle sizes ranging from 0.15 to 0.25 mm. The flat-joint model (FJM) and the smooth-joint model (SJM) are used to simulate shale matrix and layer planes, respectively. The details of these two models are exhibited by Wang et al. (2014) and Wu et al. (2018).

Calibration for the microparameters of the intact shale model
As for the PFC2D model, it is difficult to directly obtain the microparameters of the particles and the contacts through experimental results. The most common method to determine these microparameters is transformed from the macroscopic mechanical properties of samples via trial-and-error approach (Luo et al. 2018). Accordingly, the microparameters of the intact shale model are corrected (Table 1), based on the comparison between load-displacement curves, failure strength and fracture patterns.
The Brazilian test is a quasi-static test. In the numerical model, the load plates are applied on the top and the bottom of the rock specimen with infinite stiffness. The vertical displacement of the load plate is related to the loading time.
To simulate the quasi-static loading process, it is necessary to keep the load plate moving vertically at an extremely low rate during test. It should be noted that the loading rate in the numerical model is different from that in the experiment (that is, 0.05 mm/min). Basically, PFC adopts explicit time integration algorithm and its time step is determined based on the stiffness and mass of particles (Dou et al. 2019). Herein, the PFC simulation uses displacement loading with a loading rate of 0.2 m/s, which tends to be converted to 3.6 × 10 −9 m/step when the time step is set to 1.8 × 10 −8 s/ step. That is, the load plate moving 1 mm requires approximately 277,777 steps. Therefore, 0.2 m/s is slow enough to simulate quasi-static loading (Gao et al. 2016).

Comparison of the numerical and experimental results for intact shale
The testing samples were taken from fresh outcrops of the Silurian Longmaxi formation. To consider the influence of the layer angle on the failure strength, the layer angle (α) is defined as the angle between the loading direction and the normal direction of the layer plane ( Fig. 1). For the purpose of numerical model validation, seven α values (0°, 15°, 30°, 45°, 60°, 75° and 90°) are adopted in the experimental tests. The tensile strength of intact shale is estimated by the maximum load measured in the Brazilian splitting test using sample with standard size (thickness of 25 mm and diameter of 50 mm). Once the peak load is measured in the experiment, the tensile strength can be calculated as (Ma and Huang 2018;Xu et al. 2018): where F is the peak load and D and t are the diameter and thickness of the sample, respectively. It should be noted that the tensile strength calculated from Eq. (1) is appropriate for the homogeneous isotropic rocks with fracture initiating from the center of specimen. However, shale is generally characterized as heterogeneous, where the failure is not restricted to tensile damage in specimen center for most cases. Therefore, as for the Brazilian tests of shale, the value based on Eq. (1) is usually described as the failure strength, aiming to normalize the load to the diameter and thickness (Xu et al. 2018).
The load-displacement curves during Brazilian tests of the numerical and experimental performances are presented in "Appendix 1." The results show that the load-displacement curves calculated from the model are in good agreement with the experimental outputs. In addition, after  the peak is reached, both the numerical and experimental load-displacement curves drop rapidly, indicating that the shale samples have failed.
The comparison of investigations obtained from the experiments and numerical simulation methods indicates the failure strength is anisotropic and decreases with increasing α (Fig. 2). Moreover, Fig. 2 also demonstrates the failure strength from the two methods is consistent at different α values. With regard to the fracture patterns of the shale samples, the two methods also show good agreement (see Appendix 2). Accordingly, the fracture pattern of shale with low α (0° and 15°) can be regarded as a straight-line failure, where the damage mainly occurs in the matrix. As the α increases (30° and 45°), the fracture pattern changes from straight-line failure to curve-line failure. When the α increases to 60° and 75°, part of the damage develops along the layer. Finally, when the α is 90°, the damage is the pattern of straight-line failure again, in which the specimen is destroyed along the layer plane.
In conclusion, the displacement-load curve, failure strength and fracture pattern are comprehensively compared between numerical and experimental methods. Based on the good agreements above, the numerical method with the parameters of intact shale is feasible.

Evolution of microcracks in intact shale
The processes of deformation and destruction are often accompanied by the absorption, accumulation and dissipation of energy (Zhao and Xie 2008). This part mainly studies the influence of the layer angle α on the absorption energy and the damage evolution during the loading process and further reveals the mechanical characteristics of shale.
Here, the absorption energy is the energy exerted by the external force and is calculated as the area enclosed by the load-displacement curve and the abscissa axis. The absorption energy curves of intact shale with different α values are plotted in Fig. 3. In this figure, the axial displacement ratio is defined as the ratio of the loading displacement to the displacement in the peak of loading force. It can be seen that the absorption energy evolution curve grows nonlinearly. The growth rate of adsorption energy is small in the initial stage for different α values, corresponding to the compaction phase of the samples. Then, the growth rate gradually increases due to the occurrence of microcracks. When the axial displacement ratio is 80%-90%, the growth rate of energy adsorption is essentially stable (revealed by the slope of curves in Fig. 3). Besides, the distinctions among samples indicate that α has a significant influence on the absorption energy, where the sample with lower α tends to absorb more energy with a higher growth rate. Previous viewpoints suggest that the growth rate of absorption energy reflects the severity of shale failure (Zhang et al. 2018a, b), which keeps consistent with damage mechanisms of shale plotted in Fig. 4. The shale with lower α contains a dominant proportion (and a higher quantity) of microcracks in shale matrix and thus involves a greater growth rate of absorption energy, while that with higher α has an inferior proportion (and a lower quantity) of microcracks in the matrix and thus corresponds to a lower growth rate of absorption energy (Figs. 3, 4 and 5).
The damage occurring inside the shale samples can be characterized by the evolution of microcracks during the loading process, of which the entire shale failure process can be divided into three stages, as shown in Absorption energy, J ig. 3 Energy absorption curves of intact shale during loading process micropores inside the shale elastically deform and gradually close under compression. At this time, the absorption energy is converted into the elastic potential energy of the shale sample. When the displacement is loaded to point A, microcracks begin to appear (Fig. 5). During the AB stage, part of the absorption energy is converted into elastic potential energy, and the rest is released in the form of damage. Hence, the sample can still bear a load, and the rate of energy growth increases steadily (Fig. 3). In AB phase, the α has an effect on the increasing characteristics of the microcrack extent, where the development of microcracks is described in three manners. Firstly, when α is low (0° and 15°), the matrix damage determines the failure of shale samples and appears as a straight-line type. Here, the evolution characteristic shows an "incremental" mode, where the growth rate of the crack number increases gradually (Fig. 5a). Secondly, layer damage tends to occur with increasing α (30° ~ 60°), while the matrix damage still plays a dominant role. The straight-line type of failure changes to the curve-line type. In this condition, the evolution of microcracks is a "smooth" mode and rises slowly at a fixed rate (Fig. 5b). Thirdly, when the α increases to 75° and 90°, the layer damage leads to the failure of the shale samples. This evolution characteristic is described as a "stepped" mode, indicating that microcracks develop quickly in a period of time and then remain undeveloped, with repeats for several times (Fig. 5c). With regard to the BC stage, the damage develops unstably with a rapid increase of microcracks for all α (Fig. 5), in which the absorption energy and the released energy reach a relatively balanced state until the sample has completely failed. In addition, the source of microcracks (tensile or shear) in Fig. 5 keeps pace with that shown in Fig. 4, further suggesting that the fracture patterns and failure mechanisms determine the evolution of the microcracks (that is, the damage process).

Numerical model of shale embedded with fracture
According to the simulation of intact shale discussed above, the numerical approach is reliable to assess the fracture behavior of shale under Brazilian tests, in that the numerical results are in good agreement with the experimental ones. On this basis, a single fracture is embedded in the center of the specimen to study the mechanical properties of shale with natural fracture under the conditions of Brazil test. The ratio of the embedded fracture length to the sample radius is 0.3 that keeps same with some existing experimental operations and is wide adopted in numerical simulations of Brazilian tests (Haeri et al. 2014;. In this section, the influence of the relative position relationship between the fracture and layer is mainly considered, where the definition of fracture angle (β) is shown in Fig. 1.

Influence from α and β on the failure strength
It can be seen that the failure strength of shale embedded with fracture exhibits anisotropy (Fig. 6a). When β remains unchanged, the failure strength decreases along with the increasing α, which is similar to the behavior of intact shale. According to Fig. 6a, for different β values, the failure strength anisotropy of shale shows different characteristics with the α variation. When β is 0° or 15°, the influence of the embedded fracture has little effect on the failure strength and its anisotropy, in which the failure strength shows a decreasing trend over the entire α interval. As for β = 30°, the failure strength is lower than that of intact shale. Under this condition, the variation trend of failure strength decreases from α = 0° to α = 45° and then remains roughly unchanged. For higher β of 45° ~ 90°, the failure strength is further lower than that when β is 30°. Here, increasing α makes the failure strength overall decline and fluctuate within a narrow range. Besides, it can be seen from Fig. 6b that when α remains constant, the failure strength decreases from β = 0° to β = 45° and increases from β = 45° to β = 90°, exhibiting a U-shaped distribution. This phenomenon suggests that a different α-β combination triggers a variable anisotropy of failure strength.

Effect of α and β on the fracture pattern
In order to describe the fracture behavior of shale under Brazilian tests, in this study, two kinds of fracture patterns are defined: original fracture (OF) represents that the generated main crack is independent of the embedded fracture, while the induced fracture (IF) symbolizes that the development of the main crack results from the existence of the embedded fracture. The typical OF and IF types from numerical calculations are shown in Fig. 7, with the rest being tidied up in "Appendix 3." When β is 0° ~ 15°, the microcracks occur through the sample and the main crack is mainly of the OF type, which is essentially consistent with the situation of intact shale. This phenomenon indicates that the existence of embedded fracture has little effect on shale failure (β = 0° ~ 15°), in terms of the entire α interval. This fracture pattern further explains why the failure strength value is the same as that of intact shale when β is relatively small (Fig. 6). When β increases to 30°, the effect of the embedded fracture on the fracture pattern becomes obvious. In this condition, the main crack is located between the loading point and the fracture tip, which is characterized as the IF type in the α range of 0° ~ 60°. Comparatively, the fracture pattern is simultaneously comprised of OF and IF types for α of 75° or 90°. When β continues to increase (≥ 45°), the fracture pattern is regarded as IF for all α values. For the shale with an IF fracture pattern, its failure strength is lower than that of intact shale (Fig. 6). For different α and β values, the fracture patterns are summarized in Fig. 8. It can be concluded that the effect of the embedded fracture on the failure behavior tends to be noticeable when the fracture angle is high. Compared with the behavior of intact shale, shale with prefabricated fracture always produces more secondary cracks (the cracks occur at the edges of the disks) during the failure process. These secondary cracks in the developing process show a common feature that they extend toward the tip of the embedded fracture (see Appendix 3).
With regard to the damage of shale, there are four types of microcracks, namely matrix tension, matrix shear, layer tension and layer shear. As for the maximum load value, the proportion of different microcracks changes with α variation complying with a similar tendency for different β values (taking β = 15° as an example in Fig. 9). Generally, the matrix shear strength is higher than the strength of the other three types. The damage of samples during the loading process is mainly tensile (including matrix tension and layer tension), and the number of microshear cracks is very small (Fig. 9). Moreover, the microshear cracks in the matrix can even be ignored. With increasing α, layer damage is prone to occur, while the proportion of microcracks in the matrix decreases. In addition, the proportion of cracks from layer shear (or layer tension) increases under the condition of higher α value (Fig. 9).

Evolution of microcracks in shale with embedded fracture
The numerical simulations indicate the embedded fracture with a low β (0° or 15°) has a sparing effect on the evolution of microcracks in shale. For example, the microcracks in shale at β = 0° (Fig. 10a-c) develop in a way similar to that in intact shale (Fig. 5). This phenomenon results from that low β (0° or 15°) facilitates the OF failure mode (Fig. 8), where the embedded fracture seems inconsequential with the stimulated fracture. By contrast, the effect of the embedded fracture on the evolution of microcracks in noticeable for greater β (30° -90°), where the fracture pattern can be regarded as IF type for almost all the samples (Fig. 9). The representative example (β = 90°) in Fig. 10d-f suggests there are only two types of evolution trends, namely "incremental" type (e.g., α = 45°) and "stepped" type (e.g., α is 0° or 90°). In addition, Fig. 10 further indicates more shear cracks are generated at a lower β. As mentioned above, the fracture patterns and mechanisms determine the type of evolution of the microcracks.

Sensitivity analysis of layer strength parameters on shale failure
For layered shale reservoirs, the variable condition in sedimentation and diagenesis makes the mechanical properties of layer planes different (Geng et al. 2016;Feng et al. 2019). In this section, the sensitivity of layer SMJ strength on shale failure behavior is explored based on two cases in Table 2

Sensitivity of the SMJ tensile strength on shale failure (Case 1)
In Case 1, four α values and three β values are discussed in each model. The failure strength curves suggest that increasing α makes the failure strength of all samples decrease no matter what the β is, except the situation of α = 90° and β = 90° (Fig. 11). The fracture pattern and mechanism indicate that a higher layer tensile strength induces fewer secondary cracks in the model (Figs. 12, 13, 14). Moreover, for all β scope, when α is low (0° or 30°), the failure behavior of shale in all series of Case 1 is regarded as the same (Figs. 12, 13, 14), suggesting the variation in tensile strength is not sensitive to the shale failure at low α. However, with increasing α, the failure mechanisms are different among models.

Situation of β = 0°F
or β equal to 0°, the fracture pattern of each model is of OF type. As for three models, the failure mechanism is dominated by matrix tensile failure when α is 0°, and the effect of the variation in SMJ tensile strength on the failure strength is negligible at α = 0° but becomes obvious as α increases (Fig. 11a). For the LT model, when α is 30° or 60°, many layer tensile cracks occur near the loading points. Moreover, a macroscopic tensile failure plane has formed at α = 60°, which leads to a lower failure strength than that of the other two groups (Fig. 11a). When α is 90°, the failure of the LT model and the HT model is induced by the layer tension cracks and the layer shear cracks, respectively (Fig. 12). With regard to the MT model, the layer tension cracks occur near the loading points, and the layer shear cracks occur in the middle of the disk. Hence, different failure mechanisms give rise to different failure strengths. model, the failure strength of the LT model decreases, while that of the HT model remains unchanged at α = 0° (Fig. 11b). When α is 0°-60°, the matrix tension cracks play a major role in the main damage region, where more layer shear cracks occur with increasing SMJ tension strength (Fig. 13). It should be noted that when α is 90°, the failure mechanisms are obviously different for the three models. Because the SMJ tension strength is very low in the LT model, the damage mainly develops along the layer plane from the embedded fracture tip and then toward the loading point. In addition, some layer tension cracks also occur in the middle of the disk. Comparatively, with increasing SMJ tension strength, the number of layer tension cracks decreases. For the MT model, the main damage region is composed of matrix tension cracks, layer tension cracks and layer shear cracks. However, with regard to the HT model, the layer tension cracks have not occurred, while the layer shear cracks and the matrix tension cracks determine the failure mechanism and the failure strength (Fig. 13).

Situation of β = 90°B
asically, the fracture pattern is regarded as IF type for all models at β = 90° (Fig. 14). For the situation of α = 0°, the failure strengths of three models are the same, suggesting the influence of the SMJ tension strength on the failure mechanism of each model can be treated identically (Fig. 11c). However, as α increases, the failure strength of LS model and MS model declines continuously, by contrast, while that of the HT model decreases first and then increases with the minimum value at α = 60°. This variation mainly depends on the relative magnitude of the critical failure strength of different cracks. In addition, under the conditions of β = 90° and β = 45°, the SMJ tension strength has a similar influence on the damage mechanism and stimulates IF-type fractures, compared with the situation of β = 0° (OF type fractures). Based on the analysis above, the failure strength of shale is related to the tensile strength of layer planes, but its correlation degree is affected by the directions of the layer and embedded fracture (namely α and β). The failure strength can also be reflected in the fracture pattern and failure mechanism-a greater number of shear cracks correspond a higher failure strength, referring to the same failure modes.

Sensitivity of SMJ cohesion strength on shale failure (Case 2)
In Case 2, the α and β values adopted in Case 1 are also discussed in each model. As for all models in Case 2, increasing α diminishes the failure strength of shale with different β values (Fig. 15). The numerical investigations also reveal that fewer secondary cracks tend to emerge in the model with higher cohesion strength (HS model) than the one with lower cohesion strength (LS model) (Figs. 16,17,18). The details about the influence of cohesion strength variation on shale failure behavior are exhibited as below.

Situation of β = 0°A
ccording to the failure strength curves, the failure strength is the same at α = 0° and declines gradually with increasing α, among which the decrement in failure strength for LS model is obviously greater than that for HS model (Fig. 15a). For the fracture pattern, it is similar for all three models with a OF type. When α is 0°, the damage zone is located in the middle of three models and the main cracks are of the matrix tension type (Fig. 16). Nevertheless, the secondary cracks come from layer shear at lower cohesive strength (LS model), while those are basically from layer tension in HS model at α = 0°. When α increases to 30° and 60°, the lower cohesive strength (LS model) controls the damage that develops alone the layer planes, where the damage zone deviates from the middle of LS model. Comparatively, the damage zone of the MS and HS models is closer to the model center than that of LS model (Fig. 16). Moreover, fracture patterns of three models suggest a higher layer cohesive strength induces fewer layer shear cracks but more layer tension cracks, at α = 30° or 60°. When α is 90°, the failure in the LS model and HS model is induced by the layer shear cracks and the layer tension cracks, respectively. For the MS model, the layer tension cracks occur in the middle of

Situation of β = 45°I
n this situation, all the fracture patterns of all models are of IF type (Fig. 17). For all α values, there is hardly any layer tension cracking with regard to the LS model, by contrast, while more layer tension cracks tend to occur than layer shear cracks in the MS model and HS model (Fig. 17). In particular, when α is 0°, the failure mechanism is mainly matrix tension for HS model, which leads to a much higher failure strength than that of other two models (LS and MS models). When α increases to 30° and 60°, the variation in cohesion strength has limited influence on the shale failure behavior, revealed by the failure strength and fracture patterns (Figs. 15b, 17). When α is 90°, the damage develops along the layer plane from the natural fracture tip. The difference at α = 90° is that layer shear cracks occur in the LS model, whereas layer tension cracks occur in the MS model and HS model.

Situation of β= 90°T
he variation tendency of failure strength in this situation is similar to that when β is 90°; that is, increasing α is able to diminish the failure strength and the decrement in failure strength for LS model is greater than that for HS model (Fig. 15c). In terms of the fracture pattern, IF-type fracture 1 3 emerges in all models and the main damage region is located in the middle of the models (Fig. 18). Similar to the phenomenon at β = 0° or 45°, the secondary cracks mainly come from layer shear in LS model and are mainly from layer tension in MS and HS models at β = 90° (Fig. 18). From α = 0° to α = 90°, fractures originated from matrix tension occupy a decreasing proportion. Besides, for all α values, the main damage zone is in the middle of all models.

Discussion on the layer strength parameters
Based on the numerical results discussed above, the fracture pattern of the models tends to be OF type when α is 0°, suggesting the existence of embedded fracture does not affect the fracture pattern. When β increases to 45° or 90°, the influence of the embedded fracture obviously affects the fracture pattern that is regarded as IF type. In addition, if β and the layer strength remain constant, increasing α renders the layer planes more easily damaged, and the failure strength shows a decreasing trend. Besides, higher layer strength normally increases the failure strength. However, the variation in the failure strength depends on the damage mechanism. When α is 0° and β is 0° or 90°, the matrix tension cracks dominate the failure process, and the layer strength has little effect on the failure strength of the models. With regard to the failure of the layer plane, the damage can occur under tension or shear failure, which is determined by the relative magnitude of the layer tension strength and the cohesive strength. For example, when α and β are both 90°, although the damage region is in a straight-line shape, the failure mechanism is different. In contrast to the behavior of the M (MT or MS) model, a higher layer tension strength (or layer cohesion) leads to layer shear damage (or layer tension damage). Certainly, the failure strength is also different.

Conclusion
This paper establishes a numerical model in PFC2D to investigate the failure mechanism of shale containing natural fracture under Brazilian tests, which yields the following conclusions: The failure strength and fracture pattern of shale embedded with fracture are jointly determined by α and β. For different β values, the failure strength curves can be divided into three types with the variation in α. The presence of fractures reduces the anisotropy of the shale failure strength. Numerical results further indicate that the fracture pattern is dominated by IF at lower β (< 30°), while that is dominated by OF at β > 30°. The absorption energy curve of intact shale grows nonlinearly, where a lower α enables the sample to absorb more energy. Under the conditions of different α values, the microcracks in intact shale exhibit three evaluation modes, namely "incremental" mode, "smooth" mode and "stepped" mode. By comparison, as for the microcracks in shale containing natural fracture, there are only two evolution trends: the "incremental" type at α = 0° ~ 30° and the "stepped" type for higher α.
The influence from the variation in layer strength parameters on the fracture pattern is negligible at β = 0° and becomes obvious if β > 0°. If β and layer strength remain constant, the failure strength decreases with increasing α. Moreover, growing layer strength normally increases the failure strength (except when the damage mode is dominated by matrix). As for different shale models, although the failure pattern (OF or IF) may stay the same, the failure mechanism is likely to be different, thus imparting different failure strengths.