Investigations on the Directional Propagation of Hydraulic Fracture in Hard Roof of Mine: Utilizing a Set of Fractures and the Stress Disturbance of Hydraulic Fracture

Directional fracturing is fundamental to weakening the hard roof in the mine. However, due to the significant stress disturbance in the mine, principal stresses present complicated and unmeasurable. Consequently, the designed hydraulic fracture (HF) extension path is always oblique to principal stresses. Then, the HF will present deflecting propagation, which will restrict the weakness of the hard roof. In this work, we proposed an approach to drive the HF to propagate directionally in the hard roof, utilizing a set of hydraulic fractures and their stress disturbance. In this approach, directional fracturing in the hard roof is conducted via the sequential fracturing of three linear distribution slots. The disturbed stresses produced by the first fracturing (in the middle) are utilized to restrict the HF deflecting extension of the subsequent fracturing. Then, the combined hydraulic fractures constitute a roughly directional fracturing trajectory in rock, i.e., the directional fracturing. To validate the directional fracturing approach, the cohesive crack (representing rock fracture process zone (FPZ)) model coupled with the extended finite element method (XFEM) was employed to simulate the 2D hydraulic fracturing process. The benchmark of the above fracturing simulation method was firstly conducted, which presents the high consistency between simulation results and the fracturing experiments. Then, the published geological data of the hard roof in Datong coal mine (in Shanxi, China) was employed in the fracturing simulation model, with various principal stress differences (2~6MPa) and designed fracturing directions (30~60). The simulation results show that the disturbing stress of the first fracturing significantly inhibits the deflecting propagation of the subsequent fractures. More specifically, along the direction parallel to the initial minimum principal stress, the extension distance of the subsequent hydraulic fractures is 2~3 times higher than that of the deflecting HF in the first fracturing. The fracturing trajectory of the proposed direction fracturing method deviates from the designed fracturing path by only 2~14, reduced by 76%~93% compared with the traditional fracturing method utilizing a single hydraulic fracture. This newly proposed method can enhance the HF directional propagation ability more effectively and conveniently in the complex and unmeasurable stress field. Besides, this directional fracturing method can also provide references for the directional fracturing in the oil-gas and geothermal reservoir.


Introduction
Hydraulic fracturing is a fluid-solid coupled process of injecting high-pressure fluid to fracture rock and drive the fracture to propagate. Recently, hydraulic fracturing has been widely used in safe-efficient underground mining [1][2][3], stimulation of tight oil-gas reservoirs [4][5][6], geothermal exploitation in hot dry rock [3,7,8], etc. For coal mining, directional hydraulic fracturing is operative to weaken the roof structure in the mine, reducing dynamic disasters due to the hard roof sudden caving [1,2]. In detail, the hard roof is the rock layer overlying the coal seam, with prominent characteristics of high strength, large thickness, and strong integrity. As the working face advances, hard roofs are hard to cave in time and hang in the gob covering a large area. Sudden fracture and caving of large-area hanging roofs will produce shock loading, periodic roof weighting in large areas, and wind blast (large rock masses sudden caving in the mine will squeeze air to produce a storm wind). The disturbed stress due to roof bending and breaking will generate dynamic disasters [1][2][3] such as rock burst and coal-gas outburst. To avoid the above disasters, controlling the directional extension of hydraulic fracture (HF) in rock can break the hard roof along a reasonable path and transfer stress in the roadway ( Figure 1). However, owing to the strong mining-induced stresses produced by working face stoping and roadway digging, the in situ stresses present complex distribution and indeterminate directions. The designed directional propagation path of hydraulic fractures is generally oblique to the principal direction of in situ stresses [1,2]. Then, the hydraulic fracture always propagates along a deflecting trajectory, departing from the designed extension path ( Figure 1). Accordingly, the weakening of the hard roof will decrease significantly. Therefore, controlling the HF directional extension in rocks bearing complex stresses is fundamental to avoiding dynamic disasters in mine and stimulating geoenergy (such as oil-gas and geothermal energy) reservoirs.
Understanding of HF deflecting propagation mechanism is fundamental to controlling directional fracturing rocks in complex stress fields. In detail, HF deflecting propagation is very general in many situations of mining and geoenergy (oil-gas, geothermal energy, etc.) exploitation [1][2][3][4][5][6]. The HF deflecting propagation rule has been widely investigated by triaxial fracturing physical simulation tests, theoretical analyses, and numerical simulations. The previous investigations indicate that the effects on HF deflecting propagation can be mainly divided into two types. (1) In type one of HF deflecting propagation, the designed HF extension path is oblique to principal stress directions [10][11][12][13][14], changed by the strong mining-induced stresses, fracture shadowing of multistage fractures, etc. (2) Another type of effects on HF deflecting extension focuses on the anisotropic fracture resistance due to rock beddings and weak planes (such as shale beddings and coal cleats) [15,16]. The recently published experimental investigations further enrich the above impact factors, such as the initial crack's phase angle (between crack and principal stresses) and dimension, the in situ stress difference, the injection rate, and the pore pressure surrounding the fracture. More specifically, the fluctuating and moderate rise of the injection pressure is a remarkable response of HF deflecting extension [10][11][12]. With the increment of crack phase angle and injection rate, both the fracture initiation pressure and the deflecting distance improve obviously. In contrast, due to the enhancement of principal stress difference, the hydraulic fracture tends to deflect in the direction of maximum principal stress significantly [10][11][12], i.e., the fracture deflecting distance decreases. The high pore pressure zone can attract the hydraulic fracture [13,14], deflecting the fracture propagation path. For the rock involving natural weak planes, the fracture resistance along the weak plane is much lower than that of the rock matrix. Therefore, HF will deflect and propagate along the weak plane, especially in coal seams with cleats and shale containing beddings [15,16].
The above fracture deflecting extension rules support the theoretical model's development of HF deflecting propagation. The previous HF deflecting propagation models were mainly proposed by modifying the linear elastic mixedmode (I-II) fracture models, delineating the variable HF net pressure and the shear stresses due to fluid flow along with the HF. Then, the HF deflecting propagation can be quantificationally characterized roughly [17,18]. Recently, with a deeper understanding of nonlinear fracture in the rock, the nonlinear fracture models have been employed in hydraulic fracturing numerical simulations. The advances in fracturing numerical simulations mainly involve the utilization of the cohesive crack model [13,[19][20][21] and the peridynamics model [22] for representing rock fractures. Hence, some new findings of HF deflecting propagation were obtained via the above numerical simulation methods. In detail, the HF deflecting extension distance increases with the rising viscosity of fracturing fluid [13,21]. Also, the rock plastic deformations surrounding the deflecting HF enhance the fracture extension pressure and the deflecting distance [19]. Microcrack generation surrounding the hydraulic fracture tip is a primary micromechanism of HF deflecting propagation [22][23][24][25].
Based on the above HF deflecting propagation research, the directional fracturing methods, controlling HF directional extension and inhibiting HF deflecting, were proposed and divided into three categories roughly. (1) Increasing the fluid injection rate [12,21] can improve the HF deflecting radius, which will enhance the HF directional extension.
(2) Drilling dense holes or slots in rock can change local stress fields along the HF designed extension path, inducing HF directional propagation [26,27]. (3) Raising local pore pressures by continuous injection attracts the HF to propagate directionally along the designed path, restraining HF deflecting [13,14]. The above three methods provide significant references for controlling the HF directional extension in field applications but still have limitations in three aspects.
(1) Increasing the injection rate weakens HF deflecting limitedly, which is hard to overcome the control of principal stress difference on fracture deflecting. (2) The local stress fields changed by drill holes and slots are limited. For example, the stress concentration zone surrounding the drill hole is about three times the hole diameter. Therefore, directional fracturing by the above methods requires extremely dense holes or slots, complicating the engineering construction and reducing the methods' general application. (3) For low and ultralow permeability rocks, increasing the local pore pressure field by fluid injection is a long-term process, significantly reducing the engineering efficiency.
Based on the above investigations, we proposed an approach to control the HF directional propagation in the hard roof of the mine. In this approach, directional fracturing in the hard roof is conducted via the sequential fracturing of three linear distribution slots. The disturbed stresses produced by the first fracturing (in the middle) are utilized  [28][29][30], widely mentioned in the oil-gas reservoir stimulation. To date, the fracturing shadowing theories focus on the planar hydraulic fracture (i.e., mode I fracture) ( Figure 2). More specifically, the fracture shadowing theory revealed that the fluid pressure causing the extension of a hydraulic fracture also compresses the adjacent formation, mainly in the direction perpendicular to the fracture face [30] (Figure 2(b)). Then, the fluid pressure in a closed-in passive fracture will be increased in the compressed region. The magnitude of pressure increase is a function of the distance between the two fractures (passive and active fractures, whose explanation can be found in Nomenclature), the extension pressure in the active fracture, and the overlap area. Consider a two-dimensional elastic medium under the influence of two principal stresses σ 1 and σ 2 (σ 2 > σ 1 and compression positive). Suppose the medium includes a hydraulic fracture along the x-axis, with a length of 2X f (from −X f to X f ), as shown in Figure 2(c). The internal net fluid pressure is ΔP (P − σ 1 ). Under plane strain conditions, the stress distribution around the active fracture at point A (Figure 2(c)) can be computed with Equation (1) [29,30]: where σ x and σ y are the normal stress in direction-x and direction-y, τ xy the shear stresses in the plane, and r (r 1 and r 2 ) and θ (θ 1 and θ 2 ) the polar coordinates. With Equation (1), the stress variations of stress change normal to the hydraulic fracture plane surrounding the active fracture can be obtained as in Figure 2(a). For multistage fracturing in the oil-gas reservoir, the stress disturbance (i.e., fracture shadowing) of the multistage fractures will increase the fracturing pressure and the fracture deflecting extension (repulsion or attraction between fractures). The multistage fractures tend to merge due to the stress disturbance-induced deflecting propagation. Consequently, the complexity of hydraulic fractures will reduce, which will restrain the permeability increment in the oil-gas reservoir. Therefore, weakening the HF fracture shadowing (i.e., stress disturbance) is fundamental to oil-gas reservoir fracturing.
However, fracturing in underground mining is significantly different from that in oil-gas reservoirs in two main aspects. (1) Underground mining enhances the stress 3 Lithosphere disturbance and complicates the local in situ stress distribution. The surrounding rock of the mine bears complex disturbance stresses due to the roadway digging and long-wall mining. Thus, the local stress distribution of hard roofs overlying the working face is complex. In detail, the local principal stresses and their directions at different positions differ significantly from the in situ stresses in the far field. On this basis, measuring the local principal stresses of surrounding rocks in the mine remains challenging in field applications [2]. Therefore, HF deflecting propagation is general in mine, with an unpredictable extension trajectory. (2) For the HF extension space, the dimension of the hard roof in the mine is significantly small compared to the oil-gas reservoir [1][2][3][4][5][6]. Therefore, in the hard roof, due to the relatively limited fracture extension zone, the stress disturbance between hydraulic fractures will be more remarkable. Based on the above fracture shadowing theory of planar hydraulic fracture [29,30], we can infer that the deflecting hydraulic fracture also produces remarkable stress disturbances, which can be employed to conduct directional fracturing in the hard roof (introduced in the next context).

Implementation
Step. Considering the prominent characteristics of fracturing in mines, we proposed an approach to fracturing the hard roofs. In this approach, directional fracturing in the hard roof is operated by the sequential fracturing of three approximately linear distribution slots. The disturbed local stress fields of the first fracturing can induce the directional propagation of the subsequent HF ( Figure 3). To better introduce the directional fracturing method, three instructions need to be stated first. (1) The plane in Figure 3 represents the roof rock cross-section, where the HF can be simplified by a line. (2) It is difficult to determine the local principal stresses and their directions due to the prominent disturbance stresses in mine. Thus, the minimum and maximum principal stresses (σ min and σ max ) shown in Figure 3 merely represent the local stress states in rock, rather than the horizontal and vertical in situ stress in the far field. (3) The random variable θ, the angle between σ min and the designed directional HF path, represents HF propagation in a complex and unmeasurable stress field. The directional fracturing approach, utilizing a set of hydraulic fractures and their stress disturbance, is proposed to drive HF directional extension. The hard roof includes isotropic and anisotropic conditions. Rock anisotropy will complicate the application of this method, but it is not the goal in this work. The detailed steps of the proposed method are listed as follows.
Step 1. Via drilling holes and hydraulic slotting in the hard roof, three approximately linearly distributed cracks are prefabricated along the designed fracturing direction. As shown in Figure 3, σ min is the minimum principal stress, and slots A and B (B1 and B2) represent the initial cracks corresponding to the first and subsequent fracturing. The direction between the designed fracturing direction and the minimum principal stress (σ min ) is θ. As mentioned above, θ can be regarded as a variable due to the complex and unmeasurable stress fields in hard roofs. The above slots are utilized to induce the directional crack initiation of HF. Note that the stress disturbance (fracture shadowing) perpendicular to HF will decrease gradually from the fracture initiation point along the fracture propagation path [28][29][30]. Therefore, slots B1 and B2 slightly depart from the designed fracturing directions. As shown in Figure 3, d 1 < d h , whose validation will be discussed based on the simulation results in Section 5. The HF of the subsequent fracturing is supposed to initiate and propagate in a stronger disturbed stress field of the first fracturing.
For clarity, both the hydraulic slot and the initial crack mentioned in the next context represent the prefabricated cracks inducing HF directional propagation in rock. However, the above two terms are used in different contexts. Hydraulic slot and initial crack correspond to field applications and mechanical analyses, respectively.
Actively extending fracture 4 Lithosphere Step 2. The high-pressure fluid will be first injected into slot A. In the complex and unmeasurable stress fields in mine, the designed fracturing direction is always oblique to the principal stresses ( Figure 3(b)). Thus, the HF A will deflect during propagation, as shown in Figure 3(b). Though HF A presents deflecting propagation, the fracturing opening will apply compressive stresses on the rocks perpendicular to the HF, as shown in Figure 3(b). Taking the HF in Figure 3(b) as an example, the HF A generation will enhance the stress along with direction-y surrounding slots B1 and B2. Then, the stress difference (direction-x and direction-y) surrounding slots B1 and B2 will reduce. Besides, based on fracture shadowing theory [28][29][30], the compressive effect perpendicular to the fracture extension path will weaken, with the increasing distance from the fracture initiation point. Consequently, the stress disturbance of the HF A presents inhomogeneous along the extension path. When HF A propagates to a distance of about twice d 1 in direction-x, the directions of the minimum principal stress surrounding slots B1 and B2 are supposed to deflect to the designed fracturing path roughly. Therefore, the disturbed stress fields of HF A are predicted to inhibit the deflecting propagation of HF B1 and B2, i.e., improve the directional extension of the subsequent fracturing.
Step 3. The monitoring drill hole and microseismic monitoring can be employed to determine the extension distance of HF A along the direction-x. When the extension distance mentioned above reaches approximately twice d 1 , sustain the net pressure in HF A with a low fluid injection rate, which can continuously provide disturbed stresses on slots B1 and B2. Then, the high-pressure fluid is injected into slots B1 and B2 to conduct hydraulic fracturing ( Figure 3(c)). The deflecting propagation of HF B1 and B2 will be inhibited by the stress disturbance of HF, as mentioned above. Eventually, the approximately directional extension trajectory of HF A, B1, and B2 constitutes the directional fracturing path (Figure 3(c)), which is sup-posed to be roughly consistent with the designed fracturing path.
The above method is supposed to fracture the hard roof directionally in the complex and unmeasurable stress field, which will be validated by the numerical simulation of the cohesive crack model coupled with XFEM in Section 3.

Theoretical Framework of the Cohesive Crack Model Coupled with XFEM for Simulating Hydraulic Fracturing
In this work, the cohesive crack model coupled with the extended finite element method (XFEM) is employed to simulate the newly proposed directional fracturing method based on three theoretical advantages as follows. (1) HF deflecting propagation is a typical tensile-shearing mixedmode (I-II) fracturing process [18,31]. The previous investigations reveal that the development of the fracture process zone (FPZ), a microcrack zone generated at the fracture tip ( Figure 4), is a prominent fracture characteristic of rocklike materials [8,[32][33][34]. Besides, for rock-like materials, the recently published investigations [31,35] show the tensile-shearing mixed-mode (I-II) just represents the fracture bearing tensile-shearing and compressive-shearing stresses for rock-like materials. Still, at the mixed-mode (I-II) fracture tip, the deformation in FPZ and the generation of the real fracture surface are tensile. Therefore, in this work, the fracture behavior of the deflecting HF can be regarded as tensile, and the shearing stress will change the fracture extension direction. The deformation in the tensile FPZ presents a crack-like shape ( Figure 4) and has a dominant property of tension softening. To date, the cohesive crack model was widely used to delineate FPZ, over which cohesive stress, tending to close the crack, is distributed [32,36,37]. When the FPZ starts to develop, the cohesive stress decreases with the increasing crack opening 5 Lithosphere displacement, i.e., the cohesive crack presents a tensionsoftening property. Therefore, the cohesive crack can be employed to characterize the prominent features of FPZ, such as the crack-like deformation and the tensionsoftening property. Therefore, the cohesive crack model is supposed to be applicable to simulate HF deflecting propagation. (2) The extended finite element method (XFEM) [38] is an extension of the conventional finite element method based on the concept of partition of unity [39]. XFEM allows local enrichment functions to be easily incorporated into a finite element approximation, which is suitable for modeling stationary discontinuities (the limitation of the conventional finite element method), especially for the fracture extension.
(3) Though the cohesive crack model is suitable for delineating rock fracturing, for the conventional fracture propagation simulation with the cohesive element, the cohesive surfaces align with element boundaries, and the fracture propagates along with a set of predefined paths or mesh edges. In contrast, the cohesive crack model coupled with XFEM, i.e., the XFEM-based cohesive fracturing method, can be used to simulate fracture propagation along an arbitrary path in the rock-like materials since the fracture propagation is not tied to the element boundaries in a mesh. Therefore, the advantages of the cohesive crack model and XFEM can be combined in the XFEM-based cohesive fracturing method.

Modeling Approach of the Extended Finite Element
Method. The extended finite element method (XFEM) alleviates the need to create a conforming mesh, allowing easy incorporation of local enrichment functions into finite element approximation [38,39]. The special enriched functions ensure the presence of discontinuities in conjunction with additional degrees of freedom.
For fracturing simulation, the enrichment functions typically involve two parts: a shape function, the near-tip asymptotic functions that capture the singularity around the crack tip, and a discontinuous function that represents the jump in displacement across the crack surfaces. The approximation for a displacement vector function (u) with the partition of unity enrichment is expressed as in where N i ðxÞ are the usual nodal shape functions, u i the usual nodal displacement vector associated with the continuous part of the finite element solution, a i the product of the nodal enriched degree of freedom vector, HðxÞ the associated discontinuous jump function across the crack surfaces, b α i the product of the nodal enriched degree of freedom vector, and F α ðxÞ the associated elastic asymptotic crack tip functions. u i , HðxÞa i , and ∑ 4 α=1 F α ðxÞb α i are applicable to all the nodes in the model, the nodes whose shape function support is cut by the crack interior, and the nodes whose shape function support is cut by the crack tip, respectively.
The discontinuous jump function across the crack surfaces, HðxÞ, is presented as in Figure 5 and can be expressed as in where x is a sample (Gauss) point, x * is the point on the crack closest to x, and n is the unit outward normal to the crack at x * . As shown in Figure 5, the asymptotic crack tip functions are expressed as below: The cohesive-crack-like tensile-opening deformation ( where ðr, θÞ is a polar coordinate system with its origin at the crack tip. Note that θ = 0 is tangent to the crack at the tip and ffiffi r p sin ðθ/2Þ considers the discontinuity across the crack face.
3.2. Cohesive Crack Model. As mentioned above, the FPZ ahead of the tensile-shearing mixed-mode I-II fracture tip in rock can be regarded as a cohesive crack, which presents significant tensile-opening deformation and very slight shearing-slipping deformation [31,35]. Thus, we employ the traction-separation law of the cohesive crack model to simulate the hydraulic fracture propagation. Many fracturing simulation investigations call the traction-separation process of a cohesive crack as damage evolution [13,21], whereas in this work, the cohesive crack model is illustrated based on the framework of nonlinear fracture mechanics.
The cohesive crack model [41,42] indicates that the FPZ can be regarded as a cohesive crack, over which a distributed cohesive stress is acted and tends to close the cohesive crack. FPZ starts to develop when the tensile stress reaches the cohesive tensile strength, and the cohesive crack begins to open and soften at the real fracture tip (Figure 6(a)). Cohesive stress decreases with the increasing crack opening displacement (COD) (Figure 6(b)). Once the cohesive stress vanishes, COD attains the critical, the real fracture surface per unit length generates (Figure 6(b)), and FPZ develops completely.
For the traction-separation of the cohesive crack per unit length, the previous investigations [32,34,40,43] have revealed that the relationship between the decreasing cohesive stress and the increasing COD roughly follows a linear function, i.e., the linear soften relationship, as expressed in where σ t is the tensile stress applied on the cohesive crack per unit length, σ T is the cohesive tensile strength (can be regarded as the tensile strength), σ c is the cohesive stress, Δw is the crack opening displacement, and k is the slope delineating the linear softening relationship between the decreasing cohesive stress and the increasing COD. k can also be regarded as the plastic modulus, characterizing the softening-opening behavior of the cohesive crack per unit length. The expression of k is presented as in where Δw c is the critical COD. Substitution of Equation (6) into Equation (5) yields In Equation (7), Δw/Δw c represents the linear softening process, i.e., the damage evolution illustrated in some fracturing numerical simulation investigations with cohesive materials.
In this work, the XFEM-based cohesive segment method is used to simulate HF initiation and propagation [44][45][46][47]. Via introducing the phantom nodes [45][46][47], the XFEMbased cohesive segment method allows the fracture to propagate along an arbitrary path in rock-like materials rather than the element boundaries in a mesh. Phantom nodes are superposed on the original real nodes, representing the discontinuity of the crack elements ( Figure 7). For the intact element, each phantom node is completely constrained to its corresponding real node.
Phantom nodes, which are superposed on the original real nodes, are introduced to represent the discontinuity of the cracked elements, as illustrated in Figure 7. When the element is intact, each phantom node is tied to its corresponding real node. In contrast, once the element is cut through a crack, each phantom node and its corresponding real node move apart. To simulate hydraulic fracturing, the additional phantom node with pore pressure degrees of freedom is introduced on the edges of each enriched element to model the fluid flow. The fluid flow characteristics in the cracked element will be discussed in the next context.

Fluid Flow
Features in the Fracture. As mentioned above, by introducing the additional phantom node with pore pressure degrees of freedom, the tangential flow along the fracture and the normal flow across the fracture can be simulated [44]. The fluid flow patterns in the cracked elements are presented in Figure 8(a).
The tangential flow within the cracked element surfaces is modeled with a Newtonian model. The volume flow rate density vector can be expressed by where k t is the tangential permeability, i.e., the resistance to the fluid flow, ∇p is the pressure gradient along the cracked element surfaces, and d is the opening of the cracked element surfaces. Note that the tangential permeability or the resistance to flow is represented according to Reynold's equation, as in where μ is the fluid viscosity. The normal flow across the cracked element surfaces can be simulated by defining a fluid leak-off coefficient c t for the pore fluid material. The fluid leak-off coefficient delineates the pressure-flow relationship between the phantom nodes 7 Lithosphere and the cracked element surfaces. This coefficient can be regarded as the permeability of a finite layer of material on the cracked element surfaces, as shown in Figure 8(b). The fluid leak-off coefficients are expressed as below: where q t and q b correspond to the flow rates into the top and bottom surfaces of a cracked element, p i represents the pressure at the phantom node at the cracked element edge, and p t and p b are the pore pressures on the top and bottom surfaces of a cracked element.

The Constitutive Relation of Fluid-Solid Coupling in
Porous Rock-Like Materials. The rock-like material is usually a porous medium, involving the solid skeleton and the pores. For the rock whose pores are filled with fluid, i.e., the saturated rock, total stresses applied on the rock are borne by the solid skeleton and the pore fluid. The stress borne by the rock skeleton is the effective stress, which is proportionally related to the rock deformation based on Biot's theory [48]. In this work, the effective stress theory is the basis of   8 Lithosphere fluid-solid coupling simulation for hydraulic fracturing [13,44,49], as in where σ p is the effective stress tensor, σ is the total stress tensor, α is Biot's coefficient, p f is the pore pressure, and δ k is the Kronecker tensor or unit matrix. The stress equilibrium of the rock-solid phase can be expressed by the principle of virtual work, as in Equation (12) [13,44,49]: where δε is the virtual strain, t is the surface traction vector, δv is the virtual velocity vector, f is the body force vector, V is the solid volume, and S is the surface area controlled by surface traction. The continuity equation of fluid flow, representing that the fluid inflow and outflow volume is equal at one point, is expressed as in Equation (13) [12,43,44,48]: where ρ f the fluid density, ∅ is the porosity of the permeable rock, n denotes the normal outward direction of surface S, and q m is the fluid velocity vector in the porous medium.
The fluid flow behavior in porous medium follows Darcy's law: where g is the gravity acceleration vector, k is hydraulic conductivity of the porous medium, and X is a spatial coordinate vector.

Benchmark of the Cohesive Crack Model Coupled with XFEM for Simulating Hydraulic Fracturing
In this section, we selected three published hydraulic fracturing experimental cases to validate the simulation method illustrated in Section 3. More specifically, we employ the cohesive crack model coupled with XFEM to simulate the HF deflecting propagation, using the same experimental conditions and material parameters in the published works as in Table 1 and Figure 9. Then, we compare the HF propagation trajectory obtained by the numerical simulation with the experimental results to validate the numerical simulation method used in this work.
The simulation conditions and model parameters in Table 1 and Figure 9 were used to conduct hydraulic fracturing numerical simulation. The hydraulic fracture of numerical simulation was compared with the experimental results, as shown in Figure 10. Figure 10 presents that the three numerical simulation cases' deflecting hydraulic fractures are consistent with the experimental results. Therefore, the cohesive crack model coupled with XFEM for simulating hydraulic fracturing is valid.

Numerical Simulation of Directional Fracturing Utilizing a Set of Hydraulic Fractures and Their Disturbance Stresses
In this section, the numerical simulation described in Section 4 will be utilized to validate the directional fracturing method proposed in Section 2. We referred to the published geological data (including the hard roof thickness and the rock mechanical properties) [9,53] of the Datong coal mine area (Shanxi, China) to conduct numerical simulations.

Geometry Models.
The coal seam in the Datong coal mine area is approximately 240 to 350 m deep, overlaid with tight sandstone as roof rock. The roof has an average thickness of 11 m. Thus, it is difficult for these high-strength roofs to cave during initial mining due to high thickness, good integrity, and high rigidity. Based on the geological data, to characterize the local fracturing planar area of the roof, we built 2D square geometry models with a side length (L) of 11 m, as shown in Figure 11. Note that theoretical models and simulation methods are two prominent factors related to calculation correctness. These two factors have been validated via the benchmark between simulations and experiments. Therefore, the fracturing simulation is supposed to be applicable to study the HF propagation in the hard roof of practical engineering, though the size of the hard roof is much larger than the experimental specimen. Besides, the size effect is a significant research field. However, the size effect is not the goal of this work and will not be discussed in detail in this work. As this work focuses on the deflecting extension direction (i.e., the main fracture extension direction) of hydraulic fracture, the directional fracturing approach was characterized by the 2D simulation. According to the directional fracturing method mentioned in Section 2 (Figure 3), three initial cracks (slots) of 0.5 m in length were prefabricated. Then, θ were set to 30°, 45°, and 60°to simulate a random oblique intersection between the HF extension direction and the principal stress, as shown in Figure 11. To weaken the effect of the distance between slots on the simulated fracture propagation, we make the slots distributed following a unified standard. More specifically, for the numerical simulations with different θ (Figure 11) Figure 3 in Section 2, considering the disturbance of tunnelling and working face advancing for local rock in the mine, the minimum and maximum principal stresses (σ min and σ max ) applied on the 2D model ( Figure 11) represent the local stresses rather than the in situ stresses (vertical and horizontal) in the far field. Table  1: Simulation conditions and parameters of three published experimental cases.

Lithosphere
Case 1 [50] Case 2 [27] Case 3 [14] Experimental purpose Investigating HF deflecting propagation law affected by the initial fracture deviation angle and the difference of principal stresses Investigating HF directional propagation in the slotting-modified stress field Investigating HF propagation influenced by local pore pressure Experimental procedure As shown in Figure 9(a), injecting fracturing fluid to drive the prefabricated crack (A) propagate with a deviation angle of 45°A s shown in Figure 9(b), the specimen was firstly prefabricated with two slots (B), and then, the hydraulic fracturing experiment is conducted by injecting water to drive the initial crack (A) propagate As shown in Figure 9( Superscript * represents that the parameters were not supplied in the published work and were obtained via inference. Superscript # means that the parameters were not provided in the published work and were obtained by referring to similar materials in other works.

Meshing, Parameters, and Boundary Conditions.
For the property assignment of roof sandstone, the 4-node bilinear displacement and pore pressure element were employed to mesh the geometry model. Besides, the 2-node linear displacement element is utilized to simulate fracture initiation and propagation with the extended finite element method (XFEM). After meshing, the model boundaries' displacement and pore pressure conditions were set to 0 ( Figure 12). The rock mechanical and fluid (water is commonly used for roof fracturing) related parameters are listed in Table 2.

Simulation Process and Loading
Conditions. In this section, different inclination angles (θ) and the principal stresses were set for fracturing numerical simulations ( Table 3). The parameters in Table 3 were employed to simulate HF propagate in random stress fields (complex and unmeasurable). Note that the minimum principal stress was determined which refers to the in situ stress database of Chinese coal mines [53].
According to the parameters shown in Table 3, we conducted the fracturing numerical simulations with the following steps. (1) Apply and balance the principal stresses to the whole model. (2) The high-pressure fluid was injected into slot A with an injection rate of 108 L/min, driving the propagation of HF A. Since the propagation path of slot A is oblique to the minimum principal direction of in situ stresses, HF A propagates along a deflecting trajectory. When the deflecting distance of crack A along the x-direction reaches 1.5~2 times of d 1 , decrease the fluid injection rate to 0.0108 L/min. HF A will stop extending and maintain opening with net pressure and continuously produce stress disturbance on the rock surrounding HF B1 and B2.   [14,27,50]) for the three cases listed in Table 1: (a) case 1; (b) case 2; (c) case 3.

Lithosphere
high-pressure fluid with the injection rate of 108 L/min into slots B1 and B2. Then, HF B1 and B2 will extend in the disturbed stress field of HF A, directionally fracturing the hard roof in the mine.

Simulation
Results. According to the simulation method and parameters mentioned above, the simulation results can be obtained as shown in Figures 13-16. The stress contours in Figures 13-16 represent the minimum principal stress fields, where the negative and positive values represent the compressive and tensile stresses, respectively. The trajectory lines distributed in contours mean the directions of minimum principal stresses at different positions. Taking the simulation results of case B2 (in Table 3) as an example,  12 Lithosphere we will delineate the evolution of the minimum principal stress field in detail. Then, the directional fracturing method proposed in this work can be illustrated more clearly.
Step 1. After applying and balancing the principal stresses to the whole model, the minimum principal stress is distributed uniformly along the y-direction in the model, as shown in Figure 13(a).
(1) Inject high-pressure fluid into slot A and drive HF A to propagate. Due to the designed fracturing path obliquing to the principal stress, HF A begins to deflect after  Step 2-1 Step 2-2 Step 1 Stress difference evolution (directions of axis x and axis y) at point C HF B1 and B2 continuously propagate and deflect slightly compared with HF A Step 3-2: slight deflection of HF B1 and B2 Step 3-2: directional crack of HF B1 and B2 HF B1 and B2 propagate directionally in the stress field disturbed by HF A The minimum principal stress distributes uniformly along the y-axis Step 1: Principal stress equilibrium The initial crack A starts to propagate and deflect. The minimum princisal stress direction changes due to HF A propagation Step 2-1: HF A crack and starts to deflect Step 2-2: continuous propagation of deflecting HF A Pump continuously and drive HF A to propagate until the HF A exceeds the initial crack B1 and B2 in the x direction. The distribution and direction of the minimum stress field are changed remarkably      14 Lithosphere fracture initiation (Figure 13(b)). As shown in Figure 13(b), the minimum principal stresses surrounding HF A change remarkably, whereas the minimum principal stress surrounding slots B1 and B2 just deflect by 9°. (2) Continuous injection drives the deflecting HF to propagate continuously. As shown in Figure 13(b), HF A deflects entirely to the direction perpendicular to the minimum principal stress, significantly inconsistent with the designed fracturing trajectory. Figure 13(c) shows that the deflecting distance along the y-direction is only 1~1.5 m. However, the stress field in the whole model is changed dramatically due to the increasing length of fracture A. For example, the stress difference (direction-x and direction-y) at point C (located at the upper right corner of HF B2) decreases dramatically during HF A propagation, as shown in Figure 13(f). The decreased stress difference will restrain the deflecting extension of HF B1 and B2. Besides, the minimum principal stresses surrounding HF B1 and B2 change deflect by approximately 50° (  Figure 13(c)), similar to the designed inclination angle (θ = 45°shown in Figure 11(b)). Consequently, the disturbed stress field can induce the directional propagation of HF B1 and B2.
Step 3. Reduce the injection rate to maintain the net pressure in HF A. Then, the disturbed stress can continuously change the stress field surrounding HF B1 and B2. The initial cracks B1 and B2 are pumped with fracturing fluid and begin to propagate at the disturbed stress field, as shown in Figure 13(d). Since the direction of the minimum principal stress (disturbed by HF A) is almost perpendicular to the designed fracturing path, HF B1 and B2 tend to propagate along the designed path after initiation. Only when the extension distance of HF B1 and B2 along the x -direction exceeds that of fracture A, the disturbed stress of HF A decreases. Then, HF B1 and B2 start to deflect ( Figure 13(e)). But overall, the combined fracturing fractures formed by HF A, B1, and B2 can achieve the directional fracturing of the roof (Figure 13(e)).
The numerical simulation results of other cases in Table 3 are shown in (Figures 14-16). The simulated stress disturbance process in each case is similar to that shown in Figures 3-5 and will not be repeated here. At different initiation angles (θ = 30°~60°) and principal stress differences (σ max /σ min = 1:33~2), the HF deflecting propagation of the subsequent fracturing (HF B1 and B2) is significantly suppressed by the disturbed stress of the first fracturing. Consequently, the proposed directional fracturing method utilizing a set of hydraulic fractures and stress disturbance can enhance the directional fracturing ability.

Evaluation of Directional Fracturing Utilizing Combined
Fractures. As shown in Figure 17(a), the line connecting the fracture points (points where HF B1 and B2 intersect the model boundaries) can be regarded as the actual fracturing direction in the hard roof. The angle (β) between the actual fracturing direction and the designed direction can delineate the deviation of the proposed directional fracturing method.

Lithosphere
The change of β with σ max /σ min obeys different laws at different θ. In detail, when θ is 30°, β decreases with the increasing σ max /σ min (this result will be discussed in the next context). In contrast, when θ is 45°and 60°, β increases with σ max /σ min . In general, β is 2°~14°with the newly proposed method. If fracturing the roof with a single hydraulic fracture, β may reach 27°~50°(inferred with HF A). Therefore, evaluating the directional fracturing ability with β, the error of the directional fracturing method proposed in this work is just 7%~24% (i.e., reduced by 76%~93%) of that with a single hydraulic fracture.
As shown in Figure 17(b), along the direction of the initial minimum principal stress, the extension distance (w 2 ) of the HF B1 and B2 is at least 2~3 times that (w 1 ) of HF A. The higher the ratio w 2 /w 1 , the better the HF deflecting propagation can be suppressed. Note that this evaluation approach is related to the model dimension. More specifically, for the simulation case A1 in Table 3 (θ = 30°and σ max /σ min = 1:33) as in Figure 14(a2), the deflecting distance of HF A is remarkable, which significantly disturbs the HF B1 and B2. HF B1 and B2 propagate parallel to the initial minimum principal stress. This study focuses on HF propagation in the hard roof, i.e., the rock with a limited dimension. When the subsequent directional fractures (HF B1 and B2) extend to the model boundary, the value of w 2 is limited by the model size, making w 2 /w 1 in Figure 14(a2) just equal to 1. However, for the larger models, subsequent fractures HF B1 and B2 can propagate continuously and directionally, which will make w 2 /w 1 much greater than 1. Thus, case A1 in Table 3 is not presented in Figure 17(b) to avoid misunderstanding. Besides, the large β of case A1 (θ = 30°and σ max /σ min = 1:33) as shown in Figure 17(a) is also attributed to the reason mentioned above.
Above all, the proposed directional fracturing approach was preliminarily validated by the numerical simulation. This method provides references for the directional break of the hard roof in the mine. This newly proposed directional fracturing method has advantages and limitations for applications in mining. For advantages, simple construction is a remarkable advantage of the newly proposed method, including prefabricating three slots and fracturing twice in the hard roof. Besides, based on the simulation results, this method can obviously weaken the control of the complex and unmeasurable stress field on HF deflecting extension. For limitations, monitoring accuracy (by drilling monitoring holes, microseismic monitoring, etc.) of the HF extension distance during the first fracturing affects the directional extension of the subsequent fracturing. This proposed directional fracturing approach needs to be further validated in the field application. This research can also provide references for the directional fracturing in oil-gas reservoirs. For the fractured-cavity oil reservoir, oil is predominantly filled in karst caves. Consequently, directional fracturing is employed to connect karst caves and the wellbore. This directional fracturing is like "shooting a target," which can reference the directional fracturing approach proposed in this work.

Conclusions
In this work, we proposed an approach to control the HF directional propagation in the hard roof, utilizing a set of