The Impacts of Bedding Strength Parameters on the Micro-Cracking Morphology in Laminated Shale under Uniaxial Compression

Featured Application: The micro-cracking morphology and rock strength of laminated shale could be used to judge whether the shale gas formation is suitable for hydraulic fracturing or not. Abstract: The micro-cracking morphology in laminated shale formation plays a critical role in the enhancement of shale gas production, but the impacts of bedding strength parameters on micro-cracking morphology have not been well understood in laminated shale formation. This paper numerically investigated the initiation and evolution of micro-cracking morphology with bedding strength parameters in laminated shale under uniaxial compression. First, a two-dimensional particle ﬂow model (PFC 2D ) was established for laminated shale. Then, the micro-mechanical parameters of this model were calibrated using stress-strain curves and ﬁnal fracture morphology measured in the laboratory. Finally, the impacts of bedding strength parameters on the uniaxial compressive strength (UCS), crack type and the complexity of fracture network were analyzed quantitatively. Numerical simulation results indicate that the UCS of shale varies linearly with the bedding strength, especially when the shear failure of beddings is dominant. Matrix cracks mainly depend on bedding strength, while the generation of tensile cracks is determined by the shear-to-tensile strength ratio of beddings (STR). The shale with a higher STR is likely to produce a more complex fracture network. Therefore, the bedding strength parameters should be carefully evaluated when the initiation and evolution of micro-cracking morphology in laminated shale formation are simulated.


Introduction
The initiation and evolution of micro-cracking morphology in a shale gas reservoir depends not only on the hydraulic treatments but also on the mechanical properties of the rock matrix and beddings. Being different from isotropic rock, finely bedded internal structures in shale formation may result in much more complex anisotropic characteristics and thus may induce different fracturing morphologies. Previous investigations found that laminated structures could affect rock anisotropy. For instance, Cho et al. [1] tested the layered rock specimens from the area of Boryeong, Asan, Yeoncheon, and observed a significant variation of anisotropy in deformation and strength. Fjaer and Nes [2] observed dimensional numerical model within the framework of particle flow code (PFC 2D ) is built to describe the matrix-bedding microstructure of a laminated shale specimen. Particularly, different micro-constitutive laws are used to describe rock matrix and beddings. Then, the micro-mechanical parameters are calibrated using the stress-strain curves obtained from laboratory measurements. Lastly, the impacts of the tensile and shear strength of beddings on the uniaxial compression strength, crack types and crack fractal dimension of shale specimens are investigated under different dip angles. The numerical simulations reveal the role of bedding strength parameters in tensile cracking and micro-cracking morphology in shale formations. These results are helpful in understanding the relationship between bedding strength parameters and the evolution of a fracture network in shale gas formations. It also makes it possible to predict the complexity of the shale fracture network based on the final fracture morphology of shale under uniaxial compression, which can be used to optimize hydro-fracturing design siting.

Geometric Model and Loading Parameters
A uniaxial compression model for a laminated shale specimen was established. This shale model is 50 mm in diameter and 100 mm in height, which is identical to the suggestion of International Society for Rock Mechanics (ISRM) [29]. Its porosity is set as 5%, which is within the typical range of 2-10% [30,31]. Its density is 2330 kg/m 3 , which remains unchanged in simulations [32,33]. This two-dimensional model is composed of 9,223 particles (the particles are discs in PFC 2D and spheres in PFC 3D ) and approximately 20,000 contacts. These discs distribute uniformly in a range from 0.6 mm to 1.0 mm in radius. Figure 1 shows the established shale model under uniaxial compression. This model has two components: rock matrix and beddings. The rock matrix is presented in gray and the beddings are presented in black. The dip angle of beddings is defined as the angle between the bedding and horizontal plane. The bedding is horizontal if the dip angle is 0 • . In order to avoid the interference of other factors such as interlocking problems in PFC [34], the perpendicular distance of adjacent beddings maintains a constant at different dip angles of beddings. Based on previous modeling experiences [35], the perpendicular distance between adjacent beddings is set to 16 mm, which is approximately 10 times the disc diameter. The quasi-static loading rate in these numerical compression tests is set as 0.025 m/s, which refers to an investigation of loading rate effects on the cracking behavior of a rock specimen [36]. The loading process will be terminated only if the residual strength decreases to 70% of peak strength.
Appl. Sci. 2020, 10, x This study will investigate the impacts of bedding strength parameters on the micro-cracking morphology in laminated shale formations through micro-mechanical simulations. Firstly, a two-dimensional numerical model within the framework of particle flow code (PFC 2D ) is built to describe the matrix-bedding microstructure of a laminated shale specimen. Particularly, different micro-constitutive laws are used to describe rock matrix and beddings. Then, the micro-mechanical parameters are calibrated using the stress-strain curves obtained from laboratory measurements. Lastly, the impacts of the tensile and shear strength of beddings on the uniaxial compression strength, crack types and crack fractal dimension of shale specimens are investigated under different dip angles. The numerical simulations reveal the role of bedding strength parameters in tensile cracking and micro-cracking morphology in shale formations. These results are helpful in understanding the relationship between bedding strength parameters and the evolution of a fracture network in shale gas formations. It also makes it possible to predict the complexity of the shale fracture network based on the final fracture morphology of shale under uniaxial compression, which can be used to optimize hydro-fracturing design siting.

Geometric Model and Loading Parameters
A uniaxial compression model for a laminated shale specimen was established. This shale model is 50 mm in diameter and 100 mm in height, which is identical to the suggestion of International Society for Rock Mechanics (ISRM) [29]. Its porosity is set as 5%, which is within the typical range of 2-10% [30,31]. Its density is 2330 kg/m 3 , which remains unchanged in simulations [32,33]. This two-dimensional model is composed of 9,223 particles (the particles are discs in PFC 2D and spheres in PFC 3D ) and approximately 20,000 contacts. These discs distribute uniformly in a range from 0.6 mm to 1.0 mm in radius. Figure 1 shows the established shale model under uniaxial compression. This model has two components: rock matrix and beddings. The rock matrix is presented in gray and the beddings are presented in black. The dip angle of beddings is defined as the angle between the bedding and horizontal plane. The bedding is horizontal if the dip angle is 0°. In order to avoid the interference of other factors such as interlocking problems in PFC [34], the perpendicular distance of adjacent beddings maintains a constant at different dip angles of beddings. Based on previous modeling experiences35 [35], the perpendicular distance between adjacent beddings is set to 16 mm, which is approximately 10 times the disc diameter. The quasi-static loading rate in these numerical compression tests is set as 0.025 m/s, which refers to an investigation of loading rate effects on the cracking behavior of a rock specimen [36]. The loading process will be terminated only if the residual strength decreases to 70% of peak strength.

Micro-Constitutive Model for Rock Matrix
The simulation in this study is based on the interaction among adjacent discs under uniaxial compression. The motion of a disc follows the Newton's second law for the translation of mass center and the rotation around the mass center. The force-displacement law among adjacent discs is based on the contact mechanics. The contacts among the discs in a shale matrix are described through the parallel bond model (PBM), which is suitable for simulating the rock deformation and failure process [37]. Figure 2 shows the schematic diagram of this contact. Before the bond is broken, the force and moment at each contact will be updated through the following force-displacement relationship among adjacent discs: where F pb n , F pb s are the normal and shear components of contact force and M pb is the parallel-bond moment; k pb n , k pb s represent the normal and shear components of contact stiffness, and the values of these two parameters are decided by the Young moduli of contact, the sum of adjacent disc radius and the stiffness ratio listed in Table 1. A is the cross-sectional area of contact, L is the length of elastic symmetrical beam, which is equal to the sum of the radii between two adjacent discs. and I is the moment of inertia of the cross-sectional area about the neutral axis of the beam. ∆U n , ∆U s are the displacement increments in the normal and shear directions, respectively. They correspond to the extrusion/detachment and slip between adjacent discs. θ b is the relative bend-rotation increment, which corresponds to the rotation of discs.   Each contact may have two failure modes: tensile failure and shear failure. The failure is determined by the stress at the contact point. For the tensile failure, the tensile stress σ is calculated by For the shear failure, the shear stress τ is calculated by where R is the contact length in PFC 2D , and β is the moment-contribution factor. A micro-crack is generated when the contact stress exceeds the strength parameters of the bonded contact. This implies that two failure modes may occur at each contact. For example, a micro-tensile crack of rock matrix (Pbt) occurs if the tensile stress σ exceeds the tensile strength of rock matrix σ

Micro-Constitutive Model for Bedding
Bedding in shale formation is a special element. This study adopts the smooth joint model (SJM) to simulate the contacts of beddings [38]. The SJM can only deliver contact force. The tensile failure between adjacent discs in the SJM is similar to that in the PBM, but the SJM treats the shear slip between discs in a different way, as shown in Figure 3. The most distinguished feature of the SJM is that the discs along the beddings can overlap and slide across each other, while the discs can only move around one another in the PBM. Further, the dilatancy effect is considered in the shear force of the SJM if shear failure occurs. Therefore, the SJM is suitable for simulating the shear slip of beddings in rocks. The force and stress after failure in the SJM are calculated using the following fundamental algorithm:

Complexity Measure of Micro-Cracking Morphology
Micro-cracking morphology in shale formation can be comprehensively measured using the fractal dimension of the crack network. A bigger fractal dimension indicates more complex fracture morphology. Recently, fractal analysis has been successfully used to characterize the pore structures or fracture networks in an investigation of unconventional gas reservoirs [39,40]. The pore structures and fracture complexity of shale gas reservoirs control the capacity of gas flow during hydraulic fracture treatments [41]. According to the simulated final fracture morphology of the shale specimen, the effects of bedding strength parameters on fracture complexity are quantitatively analyzed through crack fractal dimension. The crack fractal dimension is calculated using the box-counting method with the following formula.
where M(n) is the number of boxes containing cracks, and 2 n represents the total number of square boxes which fully cover the binary image. The calculation procedure is stated below. A binarization algorithm is conducted first before the image segment. After that, the binary image will be covered by 2 n square boxes. By changing the size of the square box, the surface of the binary image (micro-cracking morphology after binarization treatment) can be covered by a different number of square boxes. The total number of square boxes is always an exponential of base two. Among these square boxes, there are M(n) square boxes containing cracks. As the box size decreases, the slope is the crack fractal dimension D when lgM(n) and nlg2 are projected in an orthogonal logarithmic coordinate system [42]. In this study, the crack fractal dimension is used to measure the complexity of micro-cracking morphology under different bedding strength parameters. This is to describe the complexity of the crack network.

Experiments in Uniaxial Compression
Two sets of shale specimens (0 • and 90 • ) from Longmaxi formation were physically tested in uniaxial compression (see Figure 4). Each set tested three specimens under the quasi-static loading rate of 0.05 mm/min. These specimens are 25 mm in diameter and 50 mm in height. The specimen dimensions in tests are different from those in numerical simulations. However, they have the same height/diameter ratio (2:1). Fan et al. [43] have proved that internal stress distribution is fairly even and does not change with the size of the rock specimen when the height/diameter ratio is 2:1. This implies that the dimension difference between the tested shale specimen and the simulated shale model does not affect the stress-strain responses and micro-cracking morphology.

Determination of Micro-Parameters
Some technologies such as the indentation test have been applied to measure the strength parameters of rock minerals in a local area, but it is extremely difficult to obtain the micro-parameters of shale matrix and beddings. In our simulations, all micro-parameters used in the numerical model are calibrated through a trial and error approach. The slope of the stress-strain curve, the peak strength and the macro-fracture morphology of the rock specimens were used as characteristic indicators to calibrate the strength parameters in the numerical model [44,45]. By adjusting the Young's modulus and stiffness ratio of the particle and contacts, the slope of the simulated stress-strain curve becomes parallel to the laboratory curves [46]. Further, both shear and tensile strengths at the contact are adjusted to match the UCS. The micro-parameters are finally determined when the simulation results match well with the experimental ones.
The stress-strain responses between the laboratory results (after the elimination of discreteness curve) and numerical simulations are compared in Figure 5. Figure 5a shows that the UCS, slope and the evolution process of the simulated stress-strain curves agree well with the experimental ones at the dip angle of 0°. Figure 5b shows some differences between the laboratory and numerical results when the dip angle is 90°, but the simulated peak strength still approaches the average of the experimental ones. Further, the fracture morphology of shale at 0° and 90° is similar between the experiments and simulations, as shown in Figure 6. The calibrated micro-parameters are listed in Table 1 for the rock matrix and in Table 2 for beddings. These calibrated parameters will be used in subsequent simulations. It is noted that the bedding strength is usually weaker than the shale matrix. In this specimen, the bedding strength parameters are about half of the shale matrix after calibration. With these parameters, we will explore the influences of bedding strength parameters on the stress-strain response, the percentage of different crack types and the fractal dimension of micro-cracking morphology in the next section.

Determination of Micro-Parameters
Some technologies such as the indentation test have been applied to measure the strength parameters of rock minerals in a local area, but it is extremely difficult to obtain the micro-parameters of shale matrix and beddings. In our simulations, all micro-parameters used in the numerical model are calibrated through a trial and error approach. The slope of the stress-strain curve, the peak strength and the macro-fracture morphology of the rock specimens were used as characteristic indicators to calibrate the strength parameters in the numerical model [44,45]. By adjusting the Young's modulus and stiffness ratio of the particle and contacts, the slope of the simulated stress-strain curve becomes parallel to the laboratory curves [46]. Further, both shear and tensile strengths at the contact are adjusted to match the UCS. The micro-parameters are finally determined when the simulation results match well with the experimental ones.
The stress-strain responses between the laboratory results (after the elimination of discreteness curve) and numerical simulations are compared in Figure 5. Figure 5a shows that the UCS, slope and the evolution process of the simulated stress-strain curves agree well with the experimental ones at the dip angle of 0 • . Figure 5b shows some differences between the laboratory and numerical results when the dip angle is 90 • , but the simulated peak strength still approaches the average of the experimental ones. Further, the fracture morphology of shale at 0 • and 90 • is similar between the experiments and simulations, as shown in Figure 6. The calibrated micro-parameters are listed in Table 1 for the rock matrix and in Table 2 for beddings. These calibrated parameters will be used in subsequent simulations. It is noted that the bedding strength is usually weaker than the shale matrix. In this specimen, the bedding strength parameters are about half of the shale matrix after calibration. With these parameters, we will explore the influences of bedding strength parameters Appl. Sci. 2020, 10, 5496 8 of 21 on the stress-strain response, the percentage of different crack types and the fractal dimension of micro-cracking morphology in the next section.

Micro-Parameter Value
Normal stiffness k

Micro-Cracking Behaviors of Shale under Different Bedding Parameters
The effects of dip angle and bedding strength parameters on micro-cracking morphology will be explored in the laminated shale formation. Two new indicators are introduced here: the ratio of bedding (including tensile and shear) strength to the datum value for beddings (called joint strength ratio or JSR) and the ratio of shear strength to tensile strength of beddings (called shear-to-tensile strength ratio or STR). These two indicators may comprehensively affect the micro-cracking morphology. The stress-strain and crack growth curves of shale are presented in Figure 7 when the JSR is 1.0. When the dip angle is 0 • , the axial stress increases linearly with axial strain at the pre-peak strength stage. Figure 7a shows that no micro-cracks are observed until the axial stress is at approximately 85% of the peak strength. After the peak strength of 91 MPa, the axial stress drops abruptly and shows obvious brittle failure. Simultaneously, different types of micro-cracks occur rapidly as the axial stress approaches the peak strength. This was observed by the acoustic emission of shale in the laboratory [47]. The growth rate and quantity of matrix cracks are dominant. Thus, micro-cracks occur mainly in the rock matrix. When the dip angle is 30 • , the stress-strain curve has no obvious change in shape (see Figure 7b), but the crack growth rate changes significantly. Particularly, the growth rate of shear crack in beddings increases significantly. At this time, shear cracks of beddings initiate just as the axial stress increases to approximately 60% of the peak strength. This is much earlier than that in the case of a dip angle of 0 • . When the dip angle is 60 • , the peak strength drops to 44 MPa (see Figure 7c). The shear slip of beddings initiates when the axial stress approaches about 80% of the peak strength and then increases considerably after the peak strength. In this circumstance, the damage of the rock matrix seems negligible. The stress-strain responses and the evolution of micro-cracks at the dip angle of 90 • (see Figure 7d) are similar to those at the dip angle of 0 • . The matrix failure becomes dominant again. The stress-strain curve changes little in shape but the crack growth curve is different. The shear crack of beddings occurs first and then the matrix crack appears. Therefore, the dip angle is an important parameter to affect the local stress and the evolution of micro-cracking morphology. Figure 8 summarizes the variations of UCS with the dip angle under different JSR. When the dip angle is 0 • , 30 • and 90 • , the UCS is almost the same. This implies that the dip angle has little impact. Further, the UCS is the lowest at the dip angle of 60 • but this UCS also varies with JSR. The UCS increases from 35 MPa to 44 MPa and subsequently to 53 MPa when the JSR increases from 0.8 to 1.0 and then to 1.2, respectively. This growth rate is consistent with the increase of bedding strength. This implies that the UCS of the shale sample is the most sensitive to the bedding strength when the dip angle is 60 • . Similar behaviors were observed in uniaxial compression tests [3] and Brazilian splitting tests [9]. This is because the failure pattern was changed. The micro-cracking morphology determines the UCS of shale formation. When the dip angle is 0°, 30° and 90°, the UCS is almost the same. This implies that the dip angle has little impact. Further, the UCS is the lowest at the dip angle of 60° but this UCS also varies with JSR. The UCS increases from 35 MPa to 44 MPa and subsequently to 53 MPa when the JSR increases from 0.8 to 1.0 and then to 1.2, respectively. This growth rate is consistent with the increase of bedding strength. This implies that the UCS of the shale sample is the most sensitive to the bedding strength when the dip angle is 60°. Similar behaviors were observed in uniaxial compression tests [3] and Brazilian splitting tests [9]. This is because the failure pattern was changed. The micro-cracking morphology determines the UCS of shale formation.  When the dip angle is 0°, 30° and 90°, the UCS is almost the same. This implies that the dip angle has little impact. Further, the UCS is the lowest at the dip angle of 60° but this UCS also varies with JSR. The UCS increases from 35 MPa to 44 MPa and subsequently to 53 MPa when the JSR increases from 0.8 to 1.0 and then to 1.2, respectively. This growth rate is consistent with the increase of bedding strength. This implies that the UCS of the shale sample is the most sensitive to the bedding strength when the dip angle is 60°. Similar behaviors were observed in uniaxial compression tests [3] and Brazilian splitting tests [9]. This is because the failure pattern was changed. The micro-cracking morphology determines the UCS of shale formation.

Change of Crack Types
Four crack types are identified in numerical simulations: tensile crack in matrix (Pbt), shear crack in matrix (Pbs), tensile crack in beddings (Sjt), shear crack in beddings (Sjs). This classification is based on the failure mode and location of micro-cracks. For convenient comparison, the percentage of each type of crack is the ratio of the number of those in the same crack type to the total number of micro-cracks at the end of loading process.
The influence of dip angle on the crack type is investigated. The percentage of each type is presented in Figure 9. Figure 9a shows that the percentage of shear cracks in the matrix and beddings (Pbs+Sjs) increases first and then decreases with the dip angle. As the dip angle increases from 0 • to 90 • , the percentage of shear cracks is 50.6%, 72.2%, 92.2% and 48.3%, respectively. This indicates that the percentage of each crack type varies with not only the dip angle but also the bedding strength. The first three pictures in Figure 9 show that the average ratio of the number of tensile cracks to shear cracks in the shale model is 0.628, 0.758 and 0.776 as the JSR increases from 0.8 to 1.2. This implies that the percentage of tensile cracks in the shale specimen increases with the increase of bedding strength. The impact of the JSR on the percentage of matrix cracks is presented in Figure 9d. The percentage of matrix cracks increases with the increase of bedding strength. A larger bedding strength results in a higher growth rate of the percentage of matrix cracks. For instance, when the dip angle is 60 • , the percentage of matrix cracks (Pbt+Pbs) increases from 2.5% to 6.2% and subsequently to 14.9% as the JSR increases from 0.8 to 1.0 and then to 1.2, respectively. The growth rate of matrix cracks becomes larger when the bedding strength increases in uniaxial compression loading. These simulation results indicate that the number of tensile cracks and matrix cracks increases with the increase of joint strength. During the exploration of shale gas, the increasing number of tensile cracks in the shale matrix will facilitate the shale gas transportation during hydraulic fracturing. The increase of tensile cracks (Pbt+Sjt) in fracturing behaviors could also contribute greatly to shale gas production.

Micro-Cracking Morphology
The final fracture morphology of shale is the accumulation of different micro-cracking modes. Figure 10 presents the variations of fracture morphology of shale specimens with the JSR. These figures show that both dip angle and bedding strength have significant impacts on the final fracture morphology. When the dip angle is 0 • , most of the micro-cracks occur in the rock matrix. The matrix cracks initiate from the top of the shale specimen and then cross the beddings along the loading direction. These matrix cracks coalesce into main fractures and then divert to one side of the specimen. The final fracture morphology is consistent with the laboratory observations in Figure 6a. On the micro-scale, the damage of the rock matrix results in the splitting failure of the shale specimen, while the bedding strength only has a local influence in the sub-fractures near the damaged beddings. When the dip angle is 30 • , a lot of beddings are in shear-slip failure. This limits the extension of matrix cracks from the top of the specimen. At this time, the shear cracks of beddings decrease more significantly for larger bedding strength, and the matrix cracks extend outward more evidently. This is a combined splitting and bedding sliding failure. When the dip angle is 60 • , the failure pattern of the shale specimen transforms from the combined splitting and bedding sliding failure to a pure sliding failure. Only a few matrix cracks emerge near the damaged beddings. As the JSR increases from 0.8 to 1.2, some matrix tensile cracks (Pbt) increase slightly near the sliding beddings. At the dip angle of 90 • , a combined tensile-sliding failure of beddings occurs if the JSR is 0.8. As the JSR increases from 1.0 to 1.2, the tensile crack of beddings (Sjt) decreases. The shale specimen fails gradually due to the rock matrix.
increasing number of tensile cracks in the shale matrix will facilitate the shale gas transportation during hydraulic fracturing. The increase of tensile cracks (Pbt+Sjt) in fracturing behaviors could also contribute greatly to shale gas production.  The final fracture morphology of shale is the accumulation of different micro-cracking modes. Figure 10 presents the variations of fracture morphology of shale specimens with the JSR. These figures show that both dip angle and bedding strength have significant impacts on the final fracture morphology. When the dip angle is 0°, most of the micro-cracks occur in the rock matrix. The matrix cracks initiate from the top of the shale specimen and then cross the beddings along the loading direction. These matrix cracks coalesce into main fractures and then divert to one side of the specimen. The final fracture morphology is consistent with the laboratory observations in Figure 6a. On the micro-scale, the damage of the rock matrix results in the splitting failure of the shale specimen, while the bedding strength only has a local influence in the sub-fractures near the damaged beddings. When the dip angle is 30°, a lot of beddings are in shear-slip failure. This limits the extension of matrix cracks from the top of the specimen. At this time, the shear cracks of beddings decrease more significantly for larger bedding strength, and the matrix cracks extend outward more evidently. This is a combined splitting and bedding sliding failure. When the dip angle is 60°, the failure pattern of the shale specimen transforms from the combined splitting and bedding sliding failure to a pure sliding failure. Only a few matrix cracks emerge near the damaged beddings. As the JSR increases from 0.8 to 1.2, some matrix tensile cracks (Pbt) increase slightly near The influence of the JSR on micro-cracking morphology is quantitatively analyzed by calculating the average of the crack fractal dimension at different dip angles. As the JSR increases from 0.8 to 1.2, the average fractal dimension of the matrix crack (Pbt+Pbs) is 1.207, 1.268 and 1.281 and the average fractal dimension of the tensile crack (Pbt+Sjt) is 1.225, 1.293 and 1.291, respectively. The growth rate of the fractal dimension decreases gradually with the increase of the JSR. The above analysis shows that the dip angle controls the morphology of the main fracture, while the bedding strength could change the distribution of sub-fractures. To some extent, the increase of bedding strength can enhance the complexity of the sub-fracture distribution.

Stress-Strain Responses and Micro-Cracking Evolution
The effects of shear-to-tensile strength ratio (STR) on the stress-strain responses and micro-cracking behaviors are analyzed in this section. The tensile strength of beddings σ sj c was kept to 15 MPa, while the shear strength of beddings was set to 12 MPa, 15 MPa and 18 MPa. Therefore, the STR in this section was set to 0.8, 1.0 and 1.2, respectively. Figure 11 presents the stress-strain curves and crack growth curves with dip angles of 0 • , 30 • , 60 • and 90 • when the STR is 1.0. Generally, the stress-strain curves in uniaxial compression have three stages: the linear elastic stage, the plastic deformation stage and the brittle failure stage. They are similar to those in Figure 7. This suggests that both the JSR and STR have a negligible influence on the shape of the stress-strain curves, while the dip angle of beddings makes a difference in the trend of the crack growth curves. the sliding beddings. At the dip angle of 90°, a combined tensile-sliding failure of beddings occurs if the JSR is 0.8. As the JSR increases from 1.0 to 1.2, the tensile crack of beddings (Sjt) decreases. The shale specimen fails gradually due to the rock matrix. The influence of the JSR on micro-cracking morphology is quantitatively analyzed by calculating the average of the crack fractal dimension at different dip angles. As the JSR increases from 0.8 to 1.2, the average fractal dimension of the matrix crack (Pbt+Pbs) is 1.207, 1.268 and 1.281 Appl. Sci. 2020, 10, x the sliding beddings. At the dip angle of 90°, a combined tensile-sliding failure of beddings occurs if the JSR is 0.8. As the JSR increases from 1.0 to 1.2, the tensile crack of beddings (Sjt) decreases. The shale specimen fails gradually due to the rock matrix. The influence of the JSR on micro-cracking morphology is quantitatively analyzed by calculating the average of the crack fractal dimension at different dip angles. As the JSR increases from 0.8 to 1.2, the average fractal dimension of the matrix crack (Pbt+Pbs) is 1.207, 1.268 and 1.281  stress-strain curves, while the dip angle of beddings makes a difference in the trend of the crack growth curves. Figure 12 shows that the STR alters the downward trend of the curves for the dip angle of 60°, but the UCS of the shale specimen varies little at other dip angles. When the dip angle is 60°, the UCS increases from 29 MPa to 37 MPa and subsequently to 44 MPa as the STR increases from 0.8 to 1.0 and then to 1.2, respectively. The UCS grows consistently with the increase of the STR.

Changes of Crack Type
The STR may affect the crack type. Figure 13 presents the crack types at different STRs. It can be inferred from Figure 13a that the percentage of shear cracks occurring in the matrix and beddings (Pbs+Sjs) is 51.6%, 71.8%, 98.6% and 56.4% when the dip angle is 0°, 30°, 60° and 90°, respectively. The percentage of shear cracks (Pbs+Sjs) is higher than the percentage of tensile cracks

Changes of Crack Type
The STR may affect the crack type. Figure 13 presents the crack types at different STRs. It can be inferred from Figure 13a that the percentage of shear cracks occurring in the matrix and beddings (Pbs+Sjs) is 51.6%, 71.8%, 98.6% and 56.4% when the dip angle is 0 • , 30 • , 60 • and 90 • , respectively. The percentage of shear cracks (Pbs+Sjs) is higher than the percentage of tensile cracks (Pbt+Sjt) at all dip angles. If the dip angle approaches 60 • , the number of shear cracks of the beddings becomes dominant, and the shale specimen experiences a pure bedding sliding failure. If the tensile strength of beddings (σ sj c ) is kept constant, the STR makes a difference in the percentage of each crack type. Further, the average ratio of the number of tensile cracks to shear cracks in the shale specimen increases from 0.530 to 0.608 and subsequently to 0.758 as the STR increases from 0.8 to 1.0 and then to 1.2, respectively. Therefore, the failure mode of the shale specimen is also sensitive to the STR. A larger STR can have a higher percentage of tensile failure. The growth rate of the percentage of matrix cracks is consistent with the increase of the STR. This result is a little different from that induced by the JSR.
Appl. Sci. 2020, 10, x The STR may affect the crack type. Figure 13 presents the crack types at different STRs. It can be inferred from Figure 13a that the percentage of shear cracks occurring in the matrix and beddings (Pbs+Sjs) is 51.6%, 71.8%, 98.6% and 56.4% when the dip angle is 0°, 30°, 60° and 90°, respectively. The percentage of shear cracks (Pbs+Sjs) is higher than the percentage of tensile cracks (Pbt+Sjt) at all dip angles. If the dip angle approaches 60°, the number of shear cracks of the beddings becomes dominant, and the shale specimen experiences a pure bedding sliding failure. If the tensile strength of beddings ( sj c σ ) is kept constant, the STR makes a difference in the percentage of each crack type. Further, the average ratio of the number of tensile cracks to shear cracks in the shale specimen increases from 0.530 to 0.608 and subsequently to 0.758 as the STR increases from 0.8 to 1.0 and then to 1.2, respectively. Therefore, the failure mode of the shale specimen is also sensitive to the STR. A larger STR can have a higher percentage of tensile failure. The growth rate of the percentage of matrix cracks is consistent with the increase of the STR. This result is a little different from that induced by the JSR.    Figure 14 presents the final fracture morphology of shale under different STRs. The variation of the STR hardly changes the fracture morphology if the dip angle is 0°. The failure of the shale specimen mainly occurs in the shale matrix rather than the beddings when the dip angle is 0°. When the dip angle is 30°, a higher STR promotes a cluster generation of matrix cracks rather than outward extension. The matrix cracks mainly occur at the top and the middle of the shale specimen. When the dip angle is 60°, the dominant failure is the shear slip of the beddings, and the variation of the STR only affects the number of sliding beddings. Owing to the constant tensile strength of beddings, the  The influence of the STR on the complexity of micro-cracking morphology is quantitatively analyzed through the perspective of fractal dimension. As the STR increases from 0.8 to 1.2, the average of the fractal dimension is 1.229, 1.219 and 1.268 for the matrix crack (Pbt+Pbs) and 1.186, 1.242 and 1.293 for the tensile crack (Pbt+Sjt), respectively. The STR cannot play a role in the main fracture morphology, but the increase of the STR can enhance the complexity of sub-fractures, especially when the STR is larger than 1.0.

The Difference Between Bedding Strength and Its Shear-to-Tensile Strength Ratio
Both the variation of bedding strength and the shear-to-tensile strength ratio of beddings can change the UCS, crack type and complexity of the micro-cracking morphology in shale formation. However, the difference of their contributions has not been discussed. This section will comparatively analyze their effects and identify their difference using the indicators of JSR and STR. Both the JSR and the STR are taken as 1.0 for the datum case, where the tensile strength (σ The UCS at these bedding parameters is shown in Figure 15. When the JSR increases from 1.0 to 1.2, the UCS increases linearly at the dip angle of 60 • . If the shear strength of the beddings is kept unchanged (τ sj c = 18 MPa), the variation of the STR does not affect the UCS of the shale specimen. Therefore, the bedding strength contributes more greatly to the UCS of shale. The shear strength of beddings is the essential factor affecting the UCS, especially when the shear failure of beddings becomes the dominant fracture pattern.
Appl. Sci. 2020, 10, x The influence of the STR on the complexity of micro-cracking morphology is quantitatively analyzed through the perspective of fractal dimension. As the STR increases from 0.8 to 1.2, the average of the fractal dimension is 1.229, 1.219 and 1.268 for the matrix crack (Pbt+Pbs) and 1.186, 1.242 and 1.293 for the tensile crack (Pbt+Sjt), respectively. The STR cannot play a role in the main fracture morphology, but the increase of the STR can enhance the complexity of sub-fractures, especially when the STR is larger than 1.0.

The Difference Between Bedding Strength and Its Shear-to-Tensile Strength Ratio
Both the variation of bedding strength and the shear-to-tensile strength ratio of beddings can change the UCS, crack type and complexity of the micro-cracking morphology in shale formation. However, the difference of their contributions has not been discussed. This section will comparatively analyze their effects and identify their difference using the indicators of JSR and STR.
Both the JSR and the STR are taken as 1.0 for the datum case, where the tensile strength ( The UCS at these bedding parameters is shown in Figure 15. When the JSR increases from 1.0 to 1.2, the UCS increases linearly at the dip angle of 60°. If the shear strength of the beddings is kept unchanged ( sj c τ =18 MPa), the variation of the STR does not affect the UCS of the shale specimen.
Therefore, the bedding strength contributes more greatly to the UCS of shale. The shear strength of beddings is the essential factor affecting the UCS, especially when the shear failure of beddings becomes the dominant fracture pattern.  Further, the percentages of matrix cracks and tensile cracks with different bedding strength parameters are presented in Figure 16a and Figure 16b, respectively. Compared to the datum case in  Further, the percentages of matrix cracks and tensile cracks with different bedding strength parameters are presented in Figure 16a,b, respectively. Compared to the datum case in Figure 16a For constant shear strength of beddings (τ sj c = 18 MPa), a larger STR plays an even more significant role in the increase of tensile cracks in most circumstances. Generally, the increase of bedding strength has a greater impact on the increase of the percentage of matrix cracks, while the increase of shear-to-tensile strength ratio of beddings has a greater impact on the increase of percentage of tensile cracks.
Appl. Sci. 2020, 10, x sj c τ =18 MPa) can improve the percentage of the matrix cracks (Pbt+Pbs). The JSR induces a slightly higher percentage of matrix failure than the STR. In addition, increasing the JSR (   Table 3 summarizes the average of the fractal dimension for total cracks, matrix cracks and tensile cracks at different dip angles. These results clearly demonstrate that the shale beddings with a higher STR contribute more significantly to the complexity of micro-cracking morphology.  Table 3 summarizes the average of the fractal dimension for total cracks, matrix cracks and tensile cracks at different dip angles. These results clearly demonstrate that the shale beddings with a higher STR contribute more significantly to the complexity of micro-cracking morphology.

Conclusions
In this study, a numerical model is established for the uniaxial compression of laminated shale formation. Based on this model, the impacts of the bedding strength ratio (JSR) and shear-to-tensile strength ratio (STR) of beddings on the uniaxial compressive strength (UCS), the percentage of different crack types and the crack fractal dimension of the shale specimen have been investigated at different dip angles. Particularly, the percentage of each crack type, their transformation and the micro-cracking morphology are analyzed in detail. From these simulations and analyses, the following understandings and conclusions can be drawn.
First, the numerical model established in PFC 2D can well simulate the stress-strain responses and crack growth process of shale under uniaxial compression loadings. It can successfully mimic four failure modes and their transformation in beddings and the matrix. Therefore, this numerical model can be used to simulate the mechanical characteristics and micro-cracking behaviors of layered rock materials.
Second, the bedding shear strength plays a decisive role in the UCS when the shear failure of bedding is dominant. Under this circumstance, the bedding shear strength is the essential index to change the UCS of shale. It is also found that the growth rate of UCS remains consistent with the increase in bedding strength.
Third, the dip angle plays a decisive role in the morphology of main fractures under uniaxial compression loadings, and bedding strength parameters can change the percentage of crack types and the complexity of sub-fractures in the shale formation. Therefore, the fracture morphology in shale formation depends on the local distribution of stress and strength.
Finally, micro-cracking morphology varies with the combination of tensile strength and shear strength of beddings. Their combination plays a different role in the percentage of matrix cracks and tensile cracks. A higher JSR can generate more matrix cracks. Comparatively, a higher STR may initiate more tensile cracks. The complexity of the micro-cracking morphology of shale may increase rapidly if the STR of beddings increases to some extent.