Comparative Study on Mineral-Scale Microcrack Propagation of Shale under Different Loading Methods

Key Laboratory of Shale Gas and Geoengineering, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China Innovation Academy for Earth Science, CAS, Beijing 100029, China College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing 100049, China School of Civil Engineering, Shijiazhuang Tiedao University, Shijiazhuang, Hebei 050043, China School of Mechanics and Civil Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China


Introduction
Shale gas exists in low-porosity and low-permeability shale formations. Producing a sufficient volume of multiscale crack networks is the key to enhancing the recovery of shale gas. At present, academia and industry are more concerned about the macrocrack network, and less research is focused on the propagation of microcracks. However, existing studies show that the formation of crack networks largely depends on the growth of microcracks. Minerals are the basic units of shale. Research on the physical nature of mineral-scale microcracks can reveal their propagation mechanism.
Scholars have carried out static experimental observations on the characteristics of mineral-scale microcracks after the failure of rock through the use of optical microscopes, scanning electron microscopes (SEM), computed tomography (CT), and other types of equipment. Fujii et al. [1] found that more than 90% of the tensile cracks in granite after tensile failure exist in the form of intergranular fractures, and the remaining 10% are transgranular fractures. Minerals that form microcracks at different angles to the tensile direction are unusual. Zhong et al. [2] summarized microcrack failure modes of shale in a triaxial compression experiment and categorized them into three types: tensile failure, shear failure, and slip failure. Daigle et al. [3] compared the differences between microcracks in siliceous clay and calcareous clay. Cheng and Wong [4] compared the characteristics of tensile microcracks and shear microcracks produced by the uniaxial compression of marble and analyzed the influence of mineral grains on the types of microcracks. Lan et al. [5] carried out SEM observation of microcracks before and after a triaxial compression test of shale and found that the spatial distribution of microcracks conforms to the power-law distribution, with a defined self-affine and hierarchical structure.
is study suggests that the mechanical responses of minerals have different effect on shale fracturing.
e above experiments were all static observation and analysis and did not evaluate the growth process of microcracks. SEM with an in situ loading apparatus makes it possible to observe continuous growth of microcracks. Zhao et al. [6] found in a uniaxial compression experiment of Fangshan marble that along with initiation, propagation, and connection of microcracks, the closure of microcracks also occur. Cui et al. [7,8] observed the dynamic process of microcrack initiation and propagation and discovered microscopic en echelon intermittent fractures for the first time. ey summarized the hierarchical distribution characteristics of en echelon intermittent fractures and hypothesized that the cause may be related to the heterogeneity of minerals. Zuo et al. [9] observed the continuous propagation process of main cracks and branch cracks in the microscopic region and performed finite element calculations on their fracture toughness. Dong et al. [10] found that cracks propagate in the form of microcrack connections at the microscopic scale. In other words, the microcracks near the crack tip are always in a state of competition. Microcracks connected to the crack tip become the main crack. After the main crack propagates, microcracks that are not connected to the crack tip are closed owing to stress release. Tang [11] studied the evolution of microdamage in marble under uniaxial compression and proposed a multiscale damage power-law model for cracks. Renard et al. [12] performed CT scans on the entire process of monzonite triaxial compression under 25 MPa confining pressure and used digital image correlation methods to quantitatively analyze the initiation, closure, aggregation, and penetration characteristics of microcracks at different loading stages. ese observations of the continuous growth of microcracks found many interesting phenomena that are different from the growth of macroscopic cracks. Unfortunately, none of the observations involved the distribution of minerals in the microregion, so it is impossible to clarify the influence of minerals on the growth of microcracks and thus only hypothetical explanations can be given regarding the mechanism of different observed phenomena. e in situ experimental observation of the continuous growth of microcracks based on mineral distribution is difficult. e strain and stress of the microregion also cannot be obtained easily. Scholars have developed many numerical models based on mineral distribution, such as the universal distinct element code (UDEC) particle boundary model [13,14], three-dimensional distinct element code (3DEC)-Voronoi model [15,16], distinct element model (DEM) based on three-dimensional polygons [17], particle flow code (PFC)-grain-based model (GBM) [15,18,19], and the Irazu-GBM model [20]. Li et al. [21,22] used the finite discrete element method based on grains to study factors, including boundary conditions and rock geometry, affecting the fracture growth process in granite under uniaxial compression. Lan et al. [23] carried out a uniaxial compression experiment of a mineral distribution model based on UDEC. e authors found that heterogeneity of grain size affected the distribution of tensile stress, and rocks with such grains are more likely to produce tensile microcracks than rocks with uniform grains. Li et al. [17] studied the influence of grain size on microcrack propagation based on the 3D polycrystalline discrete element method-Voronoi model. e authors reported that in specimens with larger grain sizes, the convenient path for crack propagation is along the loading direction, which results in a series of vertical tensile fractures. When the grain size was small, the microcracks were more scattered. Xu et al. [24] established a uniaxial compression experimental model based on the microstructure of coal. ey found that microcracks at the mineral interfaces occur due to heterogeneous strain of minerals. e authors concluded that the generation of tensile microcracks dominated the failure behavior of the rock under uniaxial compression. Zhou et al. [25] established PFC-GBM to study the proportion of shear and tensile cracks during crack propagation. Saadat and Taheri [26] carried out uniaxial compression experiments based on the cohesive contact model of PFC to study the formation and propagation of microcracks in the shear zone. e abovementioned mathematical simulations of microcrack propagation based on mineral distribution models were all set as uniaxial compression tests. e results suggest that tensile fractures dominated the uniaxial compression failure behavior of rock. However, hydraulic fracturing is not a simple uniaxial compression. In hydraulic fracturing, the fractures are first partially opened, and then the proppant enters to propagate the fractures. At the same time, the high pressure of the fracturing fluid pushes both sides of the fractures to move. In this process, the shale is in a complex stress state of tension, shear, and compression [27]. Observations of hydraulic fracturing cracks show that the characteristics of microcracks produced by tensile, shear, tension-shear, and compression-shear are different [28][29][30]. Other experimental studies focused on tensile and threepoint bending loading, and the observed crack network and microcrack propagation were found to be quite different from uniaxial compression [9,31].
erefore, it has important significance to study the rupture characteristics of microcracks under different loading methods. On the macroscale, some scholars have compared the peak stress, macrocrack propagation, and fracture characteristics of rocks under different loading methods [32][33][34]. However, on the mineral scale, it is essential to determine what characteristics of the microcracks are produced under different loading methods such as tension, shear, tension-shear, and 2 Advances in Civil Engineering compression-shear. It is also essential to determine the difference between microcrack growth and evolutionary form under different loading methods and what types of microcracks dominates the fracture of the rock. ese important questions need to be addressed. To reveal the influence of different loading methods on the propagation of mineral-scale microcracks, this study used the Voronoi tessellation technique to establish a cohesive zone model (CZM) based on mineral distribution. Considering the heterogeneity of minerals, various mechanical parameters were set for different minerals. Six loading methods were designed, that is, modeling tension, tension-shear, shear, and compression-shear stresses. We compared and analyzed the crack length, fractal dimension, rupture characteristics, continuous crack propagation and turning, and en echelon intermittent fracture propagation. e comparative analysis uncovered the mechanism of how different loading methods affect the growth of mineral-scale microcracks.

Voronoi Tessellation Technique.
To establish a heterogeneous shale mineral model, we added a Voronoi tessellation technique into Abaqus using a Python script. is technique contains four steps as shown in Figure 1. First, a certain number of randomly distributed control points were generated in the model sketch (Figure 1(a)). en, the adjacent control points were triangulated based on the Delaunay triangulation (Figure 1(b)). When the circumference of the constructed triangle included control points in addition to those at the vertices, the triangle was abandoned, and new triangles were formed until all control points formed vertices. ird, the circumcenters of adjacent triangles were connected (Figure 1(c)). Finally, the control points and triangles were removed and the Voronoi grains were generated (Figure 1(d)).

Numerical Calculation Principle of the Cohesion Zone
Model. CZM was proposed by Dugdale [36] and Barenblatt [37]. is model helped in solving the research difficulties caused by the crack tip stress singularity and is widely used in the study of interface delamination of composite materials [38]. Strom and Parmigiani [39] used CZM to explain the penetration and deflection behavior of cracks when passing through the bedding surface based on the stress-energy method and verified the applicability of the model.

Constitutive Equation of
Cohesive Element. e cohesive element is based on the linear elastic constitutive model of the traction-separation law. Before the cohesive element is damaged, the stress and strain satisfy the following linear elastic relationship: (1) e symbol t represents stress, ε represents strain, the subscript n represents the normal direction, and the subscripts s and t are two tangential directions perpendicular to the normal direction. For uncoupled constitutive relations, only K nn , K ss , and K tt need to be defined.

Fracture Propagation
Criteria. Initial Damage Criterion. e maximum principal stress control criterion was used in the simulation.
at is, when the normal tensile stress or tangential stress reaches the corresponding strength, it will break: where σ n is the normal stress; σ s and σ t are the shear stresses in two directions; N max is the tensile strength; and S max and T max are the shear strengths in the two directions; σ n is expressed as e cohesive element will not produce initial damage in a purely compressed state.
Damage Evolution Criterion. In this study, the damage evolution criterion was based on effective displacement, and the damage variable D was defined as where δ f m is the effective displacement at failure, taken as 0.01 in the model, δ 0 m is the effective displacement at the initial evolution of damage, δ max m , and is the maximum effective displacement during the loading process.

Principle of Acoustic Emission Simulation.
e Abaqus simulation results were extracted using a Python script, and the results were processed by MATLAB programming to simulate the acoustic emission (AE) during the loading process [40].
rough the AE location and characteristic parameters, the crack rupture condition was further described, which provided an effective means for in-depth analysis of the crack propagation process. e method to obtain the AE location can be described as follows. e node coordinates of the damaged element were extracted and the coordinates of the center point of the damaged element were calculated, which were then used as the point of the AE event location. e number of damaged elements was regarded as Advances in Civil Engineering the number of AE events. Extracting the total energy dissipated per volume of the damaged element gave the AE energy. e parameter MMIXDME, which can be used to determine the type of rupture, is defined as follows: where G n is the Type I tensile fracture energy; G s is the Type II slip fracture energy; G t is the Type III tear fracture energy.; and G T is the sum of G n , G s , and G t . e value range of MMIXDME is between 0 and 1. When MMIXDME is 0, it indicates a tensile fracture; when the value is 1, it indicates a shear fracture, and when the value is between the 0 and 1, there is a tensile-shear mixed fracture.

Model Design.
A 2D geometric model (6 mm × 6 mm) was established in Abaqus, as shown in Figure 2(a). To create a heterogeneous mineral model of shale, this 2D geometric model was divided into 722 Voronoi units using the Voronoi tessellation technique. Each Voronoi unit represented a mineral grain, and the grain size was assigned according to the study of Ji et al. [41]. Quartz, feldspar, kaolinite, and illite were randomly assigned to Voronoi units for a content of 30%, 30%, 20%, and 20%, respectively, and are represented by gray scales from high to low in Figure 2(b). Finally, the generated Voronoi units were meshed with a 0.08-mm quadrilateral grid. Zero-thickness cohesive elements were globally inserted into grain boundaries and meshes. Each mineral set contains four-node plane strain elements (CPE4) and zero-thickness cohesive elements at the boundaries of the CPE4. e mineral boundary set was created with the remaining cohesive elements at dissimilar mineral boundaries. e parameters of the CPE4 for different minerals are listed in Table 1. e parameters of cohesive elements in minerals (at boundaries of CPE4) and cohesive elements at boundaries of dissimilar minerals are shown in Table 2. It should be noted that mechanical properties of mineral boundaries depend on bond strength, and there is currently no uniform method for characterizing the bond strength of solid clay. Many experimental studies found that compared with clay minerals, mineral boundaries have a higher bulk density and stronger cementation [42][43][44][45]. erefore, the mineral boundary strength parameters were set slightly larger than those of kaolinite but smaller than those of feldspar in this model. e total number of elements in the model was 21,242, of which 6,830 were CPE4 and 14,412 were cohesive elements. We set the total loading time to 1 s, the minimum increment size to 1E-7 s, and the maximum increment size to 0.1 s.

Loading Methods.
Six loading methods were designed as shown in Figure 3, which accounted for possible stress states of hydraulic fracturing, such as tension, shear, tension-shear, and compression-shear. e mechanical boundary conditions of the model are listed in Table 3.

Results
is study successfully compared the mineral-scale microcrack path, AE, continuous crack propagation and turning, and en echelon intermittent crack propagation under different loading methods.

Crack Path.
e longer and more complex the crack produced by fracturing, the more gas-containing holes are opened, and the greater the effect on the improvement of shale gas recovery. e crack length can be measured and the complexity of the crack path can be characterized by the fractal dimension. In this study, the box dimension method was used to calculate the fractal dimension of cracks under six different loading methods. Larger fractal dimensions indicated more complicated cracks. e principle of the box dimension method is as follows. Cover the cracked area with square boxes with side length r. Some entire boxes are empty and the rest of the boxes cover a part of the crack. N(r) is the number of boxes covering the cracks; when r ⟶ 0, the fractal dimension is as follows: e geometric features of the fracture paths are shown in Figure 4. e crack length exhibited the same trend as the fractal dimension. e crack length of tension and tensionshear was longer, and the fractal dimension was larger. e (a) (b) (c) (d) Figure 1: Four steps to generate random Voronoi polygons [35].    Advances in Civil Engineering crack length of shear and compression-shear was smaller, and the fractal dimension was lower. In loading methods I-III, the crack length and fractal dimension gradually increased, whereas in loading methods IV-VI, the crack length and fractal dimension gradually decreased.

Acoustic Emission (AE).
e AE varied for crack paths under each of the six loading methods. e type of fracture, AE counts, and AE energy were extracted, as shown in Figure 5.
In contrast to the parabolic change of crack length and fractal dimension, for loading methods I-VI, the AE counts, proportion of AE counts, and proportion of AE energy of tensile fractures decreased gradually. Meanwhile, AE counts and AE energy of shear fractures increased. e AE energy of the tensile fractures first increased and then decreased in general. In loading methods I-III, the total AE counts were  6 Advances in Civil Engineering higher than those in loading methods IV-VI, but the total AE energy for I-II was lower.

Continuous Crack Propagation and Turning.
In homogeneous materials, the crack propagation path is controlled only by force; therefore, the crack grows relatively flat and smooth. However, in heterogeneous materials such as shale, the fracture path is often complicated. ese different minerals and their boundaries affect continuous propagation and turning of cracks. Figure 6 shows the final crack path (in red) under loading methods I-VI. e continuous growth and turning of cracks in different minerals and mineral boundaries are marked with arrows in different colors.

Propagation and Turning of Cracks in Kaolinite and
Illite. Among the six loading methods, the cracks were relatively flat and smooth when they propagated in the larger grains of kaolinite and illite. e length of each segment of the crack through these regions was longer, and thus these regions accounted for a higher proportion of the total crack length. When the crack tip was blocked by quartz and feldspar, the crack turned and the turning angle was low.
When small grains of kaolinite and illite were surrounded by quartz and feldspar, many small-scale inflections were produced inside the kaolinite and illite. In loading methods I-III, these inflections were relatively small. In loading methods IV-VI, there were more inflections with larger amplitudes.

Propagation and Turning of Cracks at Mineral
Boundaries. In loading methods I-III, only a small proportion of the cracks propagated along the mineral boundary. ese segments usually propagated for a short distance along the mineral boundary and then turned into a kaolinite or illite grain. e positions of the cracks along the mineral boundary in loading methods IV-VI were similar.
In loading method IV, cracks had more inflections, whereas in loading method V, crack inflections were reduced. In loading method VI, the cracks were smooth with almost no inflection, forming a shear arc.

Propagation and Turning of Cracks in Quartz and
Feldspar.
e cracks rarely extended directly through quartz and feldspar. When it occurred at the end of the crack propagation process, the crack shape was relatively straight. However, when it occurred in the middle of the crack propagation process, there were more inflections.

En Echelon Intermittent
Fractures. In addition to continuous propagation and turning, cracks were also generated, extended, and connected in a discontinuous manner, that is, en echelon intermittent cracks were created (Figure 7). e propagation process of en echelon intermittent fractures is divided into five stages.

Quartz Feldspar
Kaolinite Illite (f ) Figure 6: Continuous growth and deflection of cracks (in red) under six different loading methods. e yellow arrows refer to crack growth in large scales of kaolinite and illite, the green arrows refer to crack growth in smaller scales of kaolinite and illite surrounded by quartz and feldspar, the blue arrows show crack growth in the boundary between different minerals, and the black arrows refer to crack growth in quartz and feldspar. Stage 5: In this stage, the remaining en echelon intermittent fractures connected and the rock was broken.
e connection type of the en echelon intermittent fractures of scales under different loading methods was quite different and is summarized as follows: (1) Wing Crack Connections. After the formation of en echelon intermittent fractures, adjacent en echelon intermittent fractures were connected by generating wing cracks at both ends, and these resembled an S-shape or a reverse S-shape. e wing crack connection was evident in the connection of two largescale en echelon intermittent fractures and was also produced in the connection of small-scale en echelon intermittent fractures. Wing crack connections appeared in all six loading methods. (2) Fragment Connection. In loading methods II-VI, there were en echelon intermittent fractures that appeared to connect in the form of fragments. In loading method III, fragments were generated in the small-scale en echelon intermittent fracture connection. e fragments in loading methods II, IV, V, and VI were generated during the large-scale en echelon intermittent fracture connection. e fragment positioning for loading methods IV-VI was the same. With the increase in lateral pressure, the size of the fragments gradually decreased. e large-scale en echelon intermittent fracture in the lower part of loading method II propagated to the upper right direction, but the crack tips were not connected with other cracks. Instead, fragments were generated in the middle of the crack that connected with the other two penetrating cracks, and the en echelon intermittent fracture that extended to the upper right was abandoned by the main crack. In loading method III, the en echelon intermittent fracture in the lower position directly propagated to the upper right and did not connect with the cracks at the upper position. Consequently, the cracks at the upper position were abandoned by the main crack.

Discussion
In rock damage theory, microdamages converge to form microcracks, and then microcracks extend or connect with other microcracks to form macrocracks [46]. Based on this theory, the strain contour and AE location before the formation of macroscopic cracks in loading methods I and IV were derived as shown in Figure 8. e AE characteristics of microcracks before macrocrack formation were extracted by programming in Python, and the AE location map was drawn using MATLAB software. In tensile loading (loading method I), the deformations of different minerals were nonuniform. Owing to the lower tensile strength, kaolinite, and illite have greater strain compared to quartz and feldspar. Because of the widespread distribution of such minerals, a large number of widely distributed microcracks were generated in the rock. Most of them were tensile microcracks, whereas very few were shear microcracks. e widely distributed tensile microcracks significantly increased the potential spread range of the cracks. e cracks could extend freely or could connect to the remaining microcracks. is led to the formation of longer length and more complicated cracks. In loading method IV (shear loading), strain first occurred at both ends of the maximum shear stress line. e strain in the rest of the rock did not reach failure displacement, even when the strain of the quartz and feldspar minerals around the notch did not reach the failure displacement. ere were no widely distributed microcracks in the front of the crack propagation path; only the right-inclined and parallel tensile microcracks were generated on the right side. When the mineral at the crack tip could not be directly sheared, the crack connected with the nearest tensile microcrack. is considerably reduced the potential spread range of cracks.
Compared with loading method IV, the lateral pressure in loading methods V and VI was equivalent to directly increasing the horizontal tensile strength of the rock. However, the increment in shear strength was the normal stress multiplied by the friction coefficient that was less than 1; therefore, the increased amplitude of shear strength was much smaller than that of tensile strength. e minerals at the crack tip could be sheared directly, but the surrounding minerals had not reached its tensile strength. Lateral pressure inhibited the generation of tensile microcracks around the crack tip. As a result, the potential spread range of cracks was further reduced, and the length and complexity of the cracks were also reduced. However, the changing trend of the crack length and fractal dimension of loading methods I-III was opposite to the trend of the amount and proportion of tensile fractures. is can be explained with the following analysis. In loading method I, the direction of the maximum tensile force is horizontal. e tensile microcracks spread and tended to propagate in the vertical direction. In loading methods II and III (tensile-shear loading), the shear and tensile forces were superimposed, causing the maximum tensile force to rotate clockwise from the horizontal direction. e tensile microcracks tended to spread in the upper right direction; as a result, the main crack connected more with the tensile microcracks on the upper right of the rock.
is led to an increase in the length and complexity of the cracks. is study showed that the length and complexity of the crack cannot be simply determined by measuring the number of tensile cracks, and the spread range of the tensile microcracks is the fundamental reason.
In geology, it is believed that tortuous cracks occur due to tension, and the smooth and straight cracks occur due to shear [2].
is study demonstrated that different loading methods affecting the complexity of cracks were in the spread range of tensile microcracks. Tensile loading can not only cause bedding and fissures to fracture but can also produce tensile microcracks in widely distributed weak minerals such as kaolinite and illite. erefore, the crack has a high degree of complexity and a large spread range. Shear loading only produces parallel tensile microcracks on one side of the shear line and can only rupture the cracks on one side. e crack propagates along the direction of the maximum shear force; therefore, the tensile microcrack on one side cannot be extended easily, and the possibility of the main crack and the tensile microcrack connecting is reduced. As a result, the shear crack has a low degree of complexity and a small spread range. e differences in the continuous propagation and turning of cracks in different minerals are mainly caused by different mechanical responses of the minerals. Figure 9 shows the superposition of the AE location and crack path. Owing to the low tensile and shear strengths of weak minerals such as kaolinite and illite, microcracks easily produced tensile or shear fractures, so there was no need to connect with the surrounding tensile microcracks.  Figure 8: In loading method I (tensile loading) and loading method IV (shear loading), the strain and microfracture before the formation of macro-cracks. (a, b) Strain contour before macro-crack formation in loading methods I and IV, respectively. (c, d) In loading methods I and IV, the AE location of the microcracks distributed in the rock before the macrocrack formed. Red dots are the AE location of tensile fractures, purple dots are the AE location of shear fractures, and the colors between the two are the AE location of tensile-shear fractures. e size of the dot represents the level of AE energy. e AE energy of microfracture before the macrocrack formed is relatively small, so the AE energy magnification of this figure is larger than the others in this paper. erefore, in most cases, the crack propagation was relatively straight and the steering was relatively smooth. When weak minerals were surrounded by hard minerals such as quartz and feldspar, there were more external barriers to crack propagation, causing many types of mineral fractures. As a result, more crack inflections occurred. e tensile and shear strengths at the mineral boundary were higher than that of kaolinite and illite. erefore, in loading methods I-III, the crack only extended at the mineral boundary for a short period of time and then went back into the weak mineral. e fracture types were all tensile fractures. Loading method IV mainly produced shear fractures at the mineral boundary. When the crack deviated from the maximum shear stress line (centerline), it connected to the tensile fracture occurred in the adjacent kaolinite and illite and then returned to the shear centerline. As the lateral pressure in loading methods V and VI increased, the tensile strength of kaolinite and illite increased, and it was not easy to generate tensile fracture. erefore, the crack continued to extend along the mineral boundary in the form of a shear arc. e tensile and shear strengths of quartz and feldspar were considerably high. e tensile strength was lower than the shear strength. erefore, when cracks passed through the quartz and feldspar, they were all tensile fractures. In loading methods I and III, the situation when cracks passed through the quartz and feldspar was similar to that of breaking the last piece of strong supporting mineral, whereas in loading methods II and VI, the cracks broke the hard mineral after several inflections.
En echelon intermittent cracks were first proposed in macroscopic fractures. Cui and Han [31] observed discontinuous microcracks in experiments for the first time and explained their hierarchical structure. Figure 10 shows the AE location of small-scale en echelon intermittent crack initiation.
e small-scale en echelon intermittent cracks were all blocked by strong minerals such as quartz and feldspar, and no cracks were generated in the strong minerals. e heterogeneous strain on strong and weak minerals Kaolinite Illite (f ) Figure 9: AE location on the fracture path under six different load methods. e yellow arrows refer to the crack growth in a wide region of kaolinite and illite, the green arrows refer to the crack growth in a small region of kaolinite and illite that was surrounded by quartz and feldspar, the blue arrows refer to the crack growth in the boundary of different minerals, and the black arrows refer to the crack growth in quartz and feldspar. Red dots are the AE location of tensile fractures, purple dots are the AE location of shear fractures, and the colors between the two are the AE location of tensile-shear fractures. e size of the dot represents the level of AE energy. caused the weak minerals to reach failure displacement first, and thus, small-scale en echelon intermittent cracks were generated.
Based on this, the occurrence of large-scale en echelon intermittent cracks was also related to the barrier of strong mineral aggregates. e AE location of the second large-scale en echelon intermittent crack formation is shown in Figure 11, and there was no AE location in the strong mineral aggregates. e heterogeneous strain between the strong mineral aggregates and the weak minerals caused the weak minerals to be blocked by the strong mineral aggregates from reaching the failure displacement. Cui and Han [31] defined the position between the en echelon intermittent cracks as rock bridges. is study found that, regardless of the formation of large-scale or small-scale en echelon intermittent cracks, the rock bridges were composed of strong minerals or strong mineral aggregates.
Regardless of the scale of the en echelon intermittent crack, the spreading of cracks is related to the loading method. In loading methods I-III, en echelon intermittent cracks were mostly tensile fractures, and these were generated perpendicular to the direction of maximum tensile force. However, in loading methods IV-VI, en echelon intermittent cracks were mostly shear fractures, which were generated in the direction close to the shear centerline, and the greater the lateral pressure, the more obvious this tendency. It is worth noting that during the formation of largescale en echelon intermittent cracks in loading methods IV and V, the shear strain was concentrated at both ends of the centerline; therefore, the middle part of the rock did not reach displacement due to shear failure but reached displacement from tensile failure. As the lateral pressure increased, the length of the crack in the middle of the rock decreased. In loading method VI, no cracks occurred in the middle part of the rock, and shear fracture occurred in the upper part. Consistent with the previous analysis, lateral pressure suppressed the occurrence of tensile fracture.
Wing cracks were the most common connection types for en echelon intermittent cracks. When there were strong minerals or strong mineral aggregates in the en echelon intermittent cracks, the wing cracks bypassed the strong minerals or strong mineral aggregates to produce fragmented connections. e AE location of the fragment formation process for loading method IV is shown in Figure 12. As the shear force could not cut the strong mineral aggregates (shown in Figure 8), tensile microcracks were formed on the right side of the shear crack. e cracks connected to the tensile microcracks around the boundary of the strong mineral aggregates and eventually formed fragments. In loading methods II and III, the tensile microcracks had a large distribution range, so there were many small fragments. e lateral pressure in loading methods V and VI restrained the generation range of tensile microcracks. As lateral pressure increased, the size of the fragments decreased. e AE location of the X-shaped crack formation process for loading methods I and IV is shown in Figure 13. In loading method I, tensile microcracks were widely generated, and en echelon intermittent cracks with different directions were formed. When the en echelon intermittent cracks extended to the center, they met in the middle to form an X-shaped crack. However, in loading method IV, the en echelon intermittent crack in the middle of the rock was blocked by strong minerals and could not continue to grow. Shear microcracks that extended to the upper left formed in the middle of the crack under shear loading. e crack in the lower part propagated and connected to it, forming an X-shaped crack. e spreading range of tensile microcracks is larger than that of shear microcracks. erefore, the size of the X crack in loading method I is larger than that in loading method IV.  Following this research, in-depth studies can still be conducted in many aspects. In the future, we could set different degrees of mineral heterogeneity by ratio and size in each loading method to further confirm the role of heterogeneity and the loading method. For different loading methods, further study is needed on the similarities and differences in the en echelon intermittent crack connections when compared with connections between naturally occurring cracks, pores, and fissures. Previous studies used fractal theory to summarize the multiscale and hierarchical distribution of cracks [11], and it was found that the en echelon intermittent cracks also show the characteristics of multiscale and hierarchical distribution. e fractal study of en echelon intermittent cracks should consider the effects of loading methods, mineral heterogeneity, pores, natural fissures, and bedding. Besides, how to make the en echelon intermittent cracks more fully connected is also an important issue for the improvement of shale gas recovery. is study provides a significant reference for future research for the loading method to increase the complexity of fractures. e combination of fractal dimension and crack length per unit area demonstrated here is a promising way to evaluate the effect of fracturing in future studies.

Conclusions
In this study, the Voronoi tessellation technique was used to establish a shale CZM based on mineral distribution, and six different loading methods were applied. Based on the rock damage theory, combined with the crack path, strain contour, and AE, the initiation and propagation of microcracks under different loading methods were compared and analyzed. e following are the main conclusions drawn from this study.
(1) e essence of different loading methods affecting crack length and complexity was the spread range of tensile microcracks. Tensile loading formed a large number of widely distributed tensile microcracks in the rock, so the crack length was longer and the twists and turns were more complicated. In contrast, shear loading only produced parallel tensile microcracks on the right side of the maximum shear line, and the cracks tended to extend along the direction of the maximum shear force, resulting in a shorter crack length and low complexity. (2) e mechanical properties of various minerals and mineral boundaries are different; therefore, the mechanical responses to different loading methods are different. is leads to differences in continuous crack propagation and steering under different loading methods. (3) e formation and propagation of en echelon intermittent cracks of different scales were mainly affected by the heterogeneity of minerals and mineral aggregates. e differences in displacements at failure of different minerals led to the generation of small-scale en echelon intermittent cracks separated by strong minerals. Similarly, the differences in displacements at failure of different mineral aggregates led to the generation of large-scale en echelon intermittent cracks separated by strong mineral aggregates. (4) e spreading direction and connection form of en echelon intermittent cracks were mainly affected by the loading method. e en echelon intermittent cracks propagated perpendicular to the direction of maximum tensile stress during tensile and tensileshear loading, with various connection methods. In shear and compression-shear loading, the en echelon intermittent cracks spread near the maximum shear stress line. As the lateral pressure increased, the cracks tended to concentrate more on the centerline, and the connection method became simpler.  16 Advances in Civil Engineering

Data Availability
e data used to support this study will be made available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.