Approximately 60% of coal resources in China are broken soft coal seams, which have the characteristics of high gas content and difficulty in fracturing directly (Sun et al. 2021). The broken soft coal seams in the Hancheng block have well-developed cleats and fractures. If directly perforating and fracturing in broken soft coal seams, a large amount of fracturing fluid is easily filtered into the cleats and fractures, resulting in a fairly low contribution to the generation of hydraulic fracture (HF) (Wang et al. 2019). A short and thick fracture network is commonly formed in a local area near the wellbore, it is quite difficult to form the main fracture with a large stimulated range, and sand plugging is prone to occur. A large amount of pulverized coal is produced during hydraulic fracturing and coalbed methane (CBM) drainage. The pulverized coal easily accumulates at the pore throats along with the migration of coalbed methane, resulting in irreversible coal seam damage, thereby leading to a short stable production period of CBM wells and poor CBM exploitation benefit (Xu et al. 2018). Aiming at the problems faced by direct fracturing in broken soft coal seams, technical solutions of indirect fracturing in the roof of broken soft coal seams were proposed.

Due to the difficulty of real-time fracture monitoring of hydraulic fracturing, numerical simulation has become a powerful means to investigate the fracture propagation and evolution law of roof fracturing and optimize the fracturing plan. The numerical simulation methods of hydraulic fracturing mainly include the finite element method, discrete element method, boundary element method and meshless method. The cohesive element method is a type of finite element method that can well describe the complex fracture propagation behavior. Based on the cohesive element method, a wealth of research results on fracture propagation and evolution law have been obtained, such as layer-penetrating propagation fractures. Liu et al. established a multilayer stacking model of different coal macrolithotypes to explore fracture initiation and propagation behaviors in bright-type and dull-type coal seams (Liu et al. 2018). Guo et al. built a seepage-stress-damage coupled finite element model of sandstone‒shale layered reservoirs (Guo et al. 2017). Zhao et al. built a 3D numerical model of multilayer co-fracturing in coal measure strata and explored the impact of geological and engineering parameters on the fracture propagation law in interbedded sandstone-coal rock formations (Zhao et al. 2018). Wang et al. conducted numerical modeling and parametric studies of indirect fracturing in the roof of cataclastic or pulverized coal seams, which validated the feasibility of indirect fracturing technology to improve the CBM development effect (Wang et al. 2019). Tan et al. established a fracture vertical extension model in isotropic layered shale strata with a thick interface using the XFEM-based CZM method and investigated the influence of a thick interface on the fracture expansion behaviors (Tan et al. 2021). Considering the 3D network model of coal discontinuities based on CT data and the constitutive equation of ductile failure and mixed mode fracture, Li et al. performed numerical simulation of horizontal well indirect fracturing in the roof of crushed soft coal seams and studied the key controlling factors and their impact on the fracture network (Li et al. 2021). Chen et al. studied the geometric evolution of 3D T-shaped and offset hydraulic fractures in coal seams by using derivation of theoretical models and numerical simulation with the cohesive element method (Chen et al. 2015). The fracture morphologies involved in the previous studies included T-shaped fractures (Guo et al. 2018; Li et al. 2018; Lyu et al. 2020), inverted T-shaped fractures (Lyu et al. 2020), I-shaped fractures (Guo et al. 2017), ǂ-shaped fractures (Zou et al. 2017), and H-shaped fractures (Chen et al. 2015; Guo 2015; Sun et al. 2020). However, most of these studies focus on the initiation and propagation of one to three types of fracture morphologies. Few studies have investigated the evolution law of various fracture morphologies with different geological, rock mechanical and engineering conditions. Therefore, this issue was studied in-depth in this article.

The analysis of factors affecting the fracture parameters of coalbed methane wells is an important part of fracturing plan design. However, the influence of different factors on the fracture propagation and evolution of roof fracturing in fractured soft coal seams has not been intensively explored. This results in a lack of clear understanding of the sensitive factors and their sensitivity degree and cannot provide guidance for the design of hydraulic fracturing. Therefore, it is indispensable to perform a sensitivity analysis of hydraulic fracture parameters. The methods of influencing factor analysis mainly include the analytic hierarchy process (AHP) (Process and Saaty 2013; Yi et al. 2020), grey incidence method (Liu and Lin 2011), numerical simulation method (Zhao et al. 2018), fuzzy comprehensive evaluation method (Zuo et al. 2009; Artun and Kulga 2020), entropy weight method (Qiao et al. 2021), and multilevel comprehensive evaluation method (Hu et al. 2018; Oluwajuwon and Olugbenga 2018; Xue et al. 2021). These methods have been applied to the domain of influencing factor analysis, but there are some unavoidable problems. It needs to obtain eigenvectors and perform a consistency check in the AHP, which will lead to complex calculations and require experts to subjectively evaluate the importance of influencing factors. To pursue higher accuracy when analyzing the problem with the grey incidence method, a large amount of data is needed. When the numerical simulation method is used to analyze the main controlling factors, many assumptions are made, and the simulation results are easily distorted. It is quite difficult to calculate the subordinating degree function with the fuzzy comprehensive evaluation method. When calculating the weights of influencing factors, if only one single subjective analysis or objective analysis method is adopted, the calculated weights may deviate from the actual weights. However, relatively accurate weight values are crucial for applying the data mining method to the fracture parameter design of hydraulic fracturing. Therefore, this article combines subjective and objective analysis methods and comprehensively selects the grey relational method, entropy weight method and AHP method to study the influencing factors of fracture parameters of CBM well fracturing, aiming at conducting a simple and scientific evaluation of the factors affecting fracture parameters and obtaining a more accurate weight value.

Herein, based on the true triaxial physical simulation experiment and field test of indirect fracturing in the sandstone roof of the No. 5 broken soft coal seam in the Hancheng block, a 3D numerical model of indirect fracturing in the roof of coal seams was established by using the cohesive element method. Parameric studies were conducted. These influencing factors include the minimum horizontal principal stress difference between the roof and coal seam, coal seam permeability, injection rate and viscosity of fracturing fluid, and perforation position. The combination weight method was used to identify the dominant factors of the fracture parameters of roof fracturing. This research could provide a reference for roof fracturing plan design.

2.1 Hydraulic fracture shape analysis

Based on the field test of indirect fracturing in the sandstone roof of the No.5 coal seam of the Hancheng block, large scale true triaxial hydraulic fracturing physical simulation experiments were performed (Zhao et al. 2022). The sandstone roof layer is selected as the perforation location. The experimental results verified the feasibility of coal seam roof fracturing, as shown in Figs. 1, 2, 3 and 4. When perforating and fracturing in the sandstone roof of crushed soft coal seams, two different fracture modes occurred. One is that hydraulic fractures only extend along the interface between the sandstone roof and coal seam, forming horizontal interface fractures. The other is that the fractures penetrate through the strata interface and continue to extend downward into the coal seam. These two different fracture propagation modes are beneficial to the stimulation of crushed soft coal seams.

Photo of the test result and schematic diagram of hydraulic fractures of rock sample 1#

Photo of the test result and schematic diagram of hydraulic fractures of rock sample 2#

Photo of the test result and schematic diagram of hydraulic fractures of rock sample 3#

Photo of the test result and schematic diagram of hydraulic fractures of rock sample 4#

When the injection rate of fracturing fluid is 12 mL/min, the crack morphologies of sample 1# are mainly “干”-shaped fractures and cross fractures. The fracture morphologies are not so complex, but the effect of layer penetration is better. The fracture morphologies of sample 3# are mainly cross fractures and T-shaped fractures, and the influence of layer penetration is not obvious. When the fracturing fluid injection rate is 20 mL/min, the fracture morphology of sample 4# is dominated by “干”-shaped fractures and T-shaped fractures. The cracks of sample 2# are mainly cross fractures and T-shaped fractures, and the layer-penetrating effect is better. When the fracturing fluid viscosity is 2.5 mPa s, the crack expansion range is larger but the effect of layer penetration is not obvious. When the fracturing fluid viscosity is 15 mPa s, the crack expansion range is small, but the effect of layer penetration is obvious. Comprehensive analysis shows that the hydraulic fracture shape varies greatly under different fracturing fluid injection rates and viscosities (Fig. 5). Exploring the fracture evolution law of roof fracturing is significant for the design of fracturing plans and the optimization of on-site construction parameters, which directly affect the exploitation benefit of CBM wells.

Fracture morphologies and injection pressure curves of indirect fracturing in the sandstone roof of the No. 5 coal seam under different injection rates and viscosities

2.2 Criteria for crack initiation and propagation

Cohesive element method can not only simulate the crack initiation and propagation behavior but also simulate the normal leak-off and tangential flow of the fracturing fluid in the crack. The damage mode of the cohesive element in the ABAQUS software follows the traction–separation criterion (Li et al. 2017). The bilinear traction–separation law is the most widely used in hydraulic fracturing for its simplicity, as shown in Fig. 6. The separation is positively correlated with the traction before the damage of the cohesive element. However, once the damage initiation criterion reaches (\({\delta }^{0}\)), rock material damage can initiate following a user-defined damage evolution law, and then the traction decreases with increasing separation. When separation achieves (\({\delta }^{f}\)), the stress becomes 0, and the rock material completely fails (Li et al. 2017). The fracture initiation was judged by the quadratic stress criterion of the cohesive element (Zhu et al. 2016). The damage evolution of cohesive elements was described by the stiffness degradation (Zhao et al. 2018).

Typical traction–separation law with linear damage evolution in CZM (Confalonieri and Perego 2019; Liu et al. 2021)

2.3 Numerical modeling

Based on the field test of well A1 in the No. 5 coal seam of the Hancheng block, China, the fracture propagation model of indirect fracturing in the roof of broken soft coal seams was built, as illustrated in Fig. 7. The dimensions of the numerical model are 80 m × 40 m × 25 m. The thicknesses of the roof layer, coal seam and floor layer are 10, 5 and 10 m, respectively. The 3D roof fracturing model is embedded with a layer of vertical cohesive element and two layers of horizontal cohesive element, which are utilized to simulate the vertical fracture extension along the vertical stress direction, as well as two horizontal interface fractures between the coal seam and the top and bottom sandstone layers. The element type of hydraulic fracture is the COH3D8P element, which can simulate the tangential flow and normal leak-off of fracturing fluid during hydraulic fracturing. The element type of the rock mass is the C3D8P element, which can simulate the 3D deformation and internal fluid seepage process of the rock mass (Zhao et al. 2018). The viscosity and injection rate of the fracturing fluid are 15 mPa s and 3 m³/min, respectively. The injection point of fracturing fluid in the numerical model is 2.5 m away from the upper interface of the coal seam.

3D numerical model of roof fracturing

Uniaxial compression tests, triaxial compression tests and fracture toughness tests were performed on natural coal and sandstone core plugs. Figure 8 shows the fracture morphology comparison of core samples before and after the tests. The test results are illustrated in Table 1.

Photos of rock samples before and after rock mechanical experiments

According to the test results of rock mechanical properties in the laboratory (Table 1) and field fracturing construction data, the modeling parameters are determined, as listed in Table 2.

2.4 Fracture initiation and extension process during roof fracturing

To observe the extension and evolution behaviors of cracks during hydraulic fracturing more directly, the fracturing time is normalized. The normalized fracturing time refers to the ratio of the real time at a certain moment to the total fracturing time. Figure 9 shows the fracture morphologies at different fracturing times. At the beginning of the fracturing simulation, local pressure accumulation was formed, resulting in pressure suppression. The continuous injection of fracturing fluid promotes vertical downward expansion of the fracture. When reaching the interface between the sandstone roof and coal seam, because of the difference in rock mechanical properties between layers, the pressure holding phenomenon was formed at the interface again. When the hydraulic pressure reaches a critical pressure, the hydraulic fracture will break through the lithological interface and propagate into the coal seam. Subsequently, hydraulic fractures mainly expand in the coal seams and have difficulty breaking through the lower sandstone floor. Finally, a “士”-shaped fracture dominated by horizontal fractures was generated. The maximum length, width and height of the horizontal fracture in the upper lithological interface are 30.06 m, 3.52 mm and 10.18 m, respectively. The maximum length, width and height of the vertical fracture in the coal seam are 25 m, 4.19 mm and 5 m, respectively. A local amplification profile of the crossing of the “士”-shaped fracture is exhibited in Fig. 10.

Fracture extension and evolution during indirect fracturing

Local magnification at the intersection of the “士”-shaped fracture

2.5 Verification by field test

The pressure data collected during the on-site fracturing operation are the pump pressure of the ground injection pump, while the pressure data calculated by the numerical simulation are the injecting point pressure underground. To compare the pump pressure curve in the numerical simulation with that of field fracturing construction, the injection point pressure in the numerical simulation can be converted into the pump pressure of the ground injection pump by Eq. (1) (Zhu et al. 2013).

where

Figure 11 shows that the on-site construction pressure declines sharply after increasing dramatically to the peak in the initial stage. During this time, the sandstone layer breaks up, and the hydraulic fracture extends in the vertical direction. After the pressure decreases to 13.55 MPa, it increases rapidly to the maximum pump pressure of 31.34 MPa and then decreases sharply, indicating that the hydraulic fracture penetrates the lithologic interface between the sandstone roof and coal seam and then continues to expand downward into the coal seam. Then the on-site pump pressure fluctuates slightly but is relatively stable as a whole, suggesting fracture propagation in the crushed soft coal seam. From the on-site pump pressure curve and the injection pressure curve of the numerical simulation, it can be inferred that the hydraulic fracture mainly expands in the coal seams during the fracturing period, and the indirect fracturing effect of crushed soft coal seams is desirable.

Pressure curves of the on-site construction and numerical simulation of well A1

When the hydraulic fracture extends in the coal seam, the change trends of the on-site pressure curve and the injection pressure curve of the numerical simulation are basically similar. The average pressure difference between the field and the numerical simulation is 2.68 MPa, and the average error is 17.98%. Adding sand was not considered in the 3D model of roof fracturing. The actual formation distribution in the field is a superposition combination in the spatial distribution. However, it is assumed that the sandstone roof and floor and coal seam are homogeneous lithologies in the numerical simulation. Hence, the error of the numerical simulation results is within a reasonable range. In conclusion, the numerical model of roof fracturing in broken soft coal seams is reliable.

3.1 Fracture propagation and evolution behaviors under different influencing factors

3.1.1 Minimum horizontal principal stress difference

The minimum horizontal principal stress difference between the sandstone roof and coal seam is a vital factor controlling the fracture vertical extension in the roof (Wang et al. 2019). The minimum horizontal principal stress of the broken soft coal seam is 6.9 MPa. By changing the minimum horizontal principal stress of the sandstone roof layer, the stress difference between layers (equal to the minimum horizontal principal stress of the sandstone layer minus that of the coal seam, denoted as \({\Delta \sigma }_{rc}\)) is −1, 0, 1, 2 and 3 MPa, respectively.

Figure 12 shows that different crack shapes are formed under different stress combinations. With the increase in the minimum horizontal principal stress difference between layers from −1 to 1 MPa, the hydraulic fracture morphology transforms from a T-shaped fracture to a cross-shaped fracture. When \({\Delta \sigma }_{rc}\) increases from 1 to 2 MPa, the hydraulic fracture changes from a cross-shaped fracture to a “士”-shaped fracture. Finally, a “士”-shaped fracture dominated by horizontal fracture is formed. Because cracks tend to extend in the direction of less resistance, the increase in the minimum horizontal principal stress difference between layers can promote the vertical downward propagation of hydraulic fractures and limit the crack height in the roof and coal seam.

Fracture morphologies under different stress differences between layers

Figure 13 shows the fracture shape parameters in three directions. Figure 14a, b illustrate that when the \({\Delta \sigma }_{rc}\) difference between layers increases from −1 MPa to 2 MPa, the height and width of the vertical fracture in the coal seam gradually increase, forming a well-developed “士”-shaped fracture. When \({\Delta \sigma }_{rc}\) is 2 MPa, the vertical fracture has already penetrated the whole coal seam, and the lower horizontal interface fracture startes to expand. Hence, when \({\Delta \sigma }_{rc}\) is higher than 2 MPa, the changes in the height and width of vertical fractures in the coal seam are not obvious (Fig. 14a), and the half-length of vertical fractures in the coal seam has a decreasing trend (Fig. 14b). Figure 14c, d show that when \({\Delta \sigma }_{rc}\) changes from 2 to 3 MPa, the width of the upper horizontal fracture has a decreasing trend. Therefore, the extension effect of hydraulic fractures is considered to be the most sufficient when \({\Delta \sigma }_{rc}\) is equal to 2 MPa. The stress difference between layers of 2 MPa can trigger desirable hydraulic fractures.

Fracture parameter labeling map in different directions

Hydraulic fracture profiles under different minimum horizontal principal stress differences between layers

3.1.2 Coal seam permeability

Coal seam permeability is one of the key influencing factors on fracture extension in coal seams during fracturing operations (Zhao et al. 2019). Therefore, five numerical simulation cases were conducted under coal seam permeabilities of 5, 9, 15, 20 and 30 mD, respectively. The simulated fracture shapes are exhibited in Fig. 15.

Fracture morphologies under different coal seam permeabilities

Figures 15 and 16 show that when the coal seam permeability increases from 5 to 15 mD, the hydraulic fracture morphology changes from a “士”-shaped fracture to a cross-shaped fracture, the height and width of vertical fractures in the coal seam gradually decrease, and the expansion range of hydraulic fractures in the coal seam gradually decreases. When the permeability of the coal seam increases from 15 to 30 mD, the hydraulic fracture morphology changes from a cross-shaped fracture to a T-shaped fracture, and the width and height of vertical fractures in the coal seam decrease significantly. The leak-off of the fracturing fluid at the upper horizontal coal seam interface is serious, resulting in the maximum height, width and length of the hydraulic fracture at the upper horizontal interface, which cannot penetrate the interface to vertically extend. It is concluded that when the coal seam permeability is higher than 15 mD, the leak-off of fracturing fluid into the upper horizontal interface of the coal seam increases, rendering its slight contribution to the vertical fracture extension, thus resulting in insufficient communication of hydraulic fractures in the coal seam.

Hydraulic fracture profiles under different coal seam permeabilities

3.1.3 Fracturing fluid injection rate

Fracturing fluid injection rate greatly affects fracture expansion during the on-site fracturing construction process (Guo et al. 2017; Wang et al. 2019). When the fracturing fluid injection rate is 1 m³/min, 2 m³/min, 3 m³/min, 4 m³/min and 5 m³/min, the simulated results are shown in Fig. 17.

Fracture morphologies under different fracturing fluid injection rates

Figures 17 and 18 show that when the injection rate of the fracturing fluid changes from 1 to 3 m³/min, the hydraulic energy is sufficient to penetrate the lithological interface and expand downward into the coal seam. Consequently, the hydraulic fracture shape transforms from a T-shaped fracture to a cross-shaped fracture. With a further increase in the injection rate, the hydraulic energy is enough to open the lithologic interface between the coal seam and floor, forming a lower horizontal interface fracture. Therefore, the hydraulic fracture morphology changes into a “士”-shaped fracture. The width, length and height of the vertical fracture in the coal seam increase significantly, and the communication range of the “士”-shaped fracture in the coal seam increases gradually. Note that when the injection rate of the fracturing fluid is 1 m³/min, the hydraulic fracture initiates from the roof and then extends along the lithological interface between the roof and coal seam. Compared with the hydraulic fractures of 2 m³/min, 3 m³/min, 4 m³/min and 5 m³/min, the hydraulic fractures of 1 m³/min did not penetrate the lithological interface to expand into the coal seam but formed the shortest vertical fracture (Fig. 17). The hydraulic energy was mainly used to create the horizontal interface fracture. In addition, the length and height of the horizontal fracture were the smallest (Fig. 18c, d). Therefore, there is sufficient hydraulic energy in the fracture width direction, resulting in the largest horizontal crack width. Figure 18a, d show that when the fracturing fluid injection rate changes from 3 to 5 m³/min, the height and length of the upper horizontal interface fracture of the “士”-shaped fracture increase significantly, and the increase in fracturing fluid flow rate is mainly used to create the horizontal interface fracture in the upper coal seam. Figure 18a, b indicate that when the injection rate of the fracturing fluid increases from 3 to 5 m³/min, the width and length of the vertical fracture in the coal seam have a decreasing trend. Therefore, during sandstone roof fracturing operations, determining the critical fracturing fluid injection rate is essential to optimize the fracturing construction scheme.

Hydraulic fracture profile morphologies under different fracturing fluid injection rates

3.1.4 Fracturing fluid viscosity

Fracturing fluid viscosity has a great impact on the leak-off of fracturing fluid and fracture extension during hydraulic fracturing (Zhao et al. 2018). Therefore, numerical simulations were conducted under fracturing fluid viscosities of 5, 10, 15, 20 and 30 mPa s, respectively. The fracture morphologies and profiles are exhibited in Figs. 19 and 20.

Fracture morphologies under different fracturing fluid viscosities

Hydraulic fracture profiles under different fracturing fluid viscosities

According to Figs. 19 and 20, when the fracturing fluid viscosity changes from 5 mPa s to 15 mPa s, the hydraulic fracture morphology changes from a cross-shaped fracture to a “士”-shaped fracture, the width and height of the vertical fracture in the coal seam gradually increase, and the height of the upper horizontal fracture of the coal seam increases significantly. When the fracturing fluid viscosity increases from 15 mPa s to 30 mPa s, the hydraulic fracture morphology transforms from a “士”-shaped fracture to a cross-shaped fracture, the vertical fracture height and width of the coal seam decrease significantly, and the horizontal fracture height and length of the upper coal seam also have a decreasing trend. When the fracturing fluid viscosity is 15 mPa s, the propagation range of the “士”-shaped fracture in the coal seam reaches the maximum. When the viscosity of the fracturing fluid is small, with the increase in the fracturing fluid viscosity, the leak-off in the fracturing process is gradually weakened, which has a promoting effect on the vertical extension of hydraulic fractures, and the communication range of the hydraulic fracture in the coal seam is gradually increased. After the fracturing fluid viscosity reaches a critical value of 15 mPa s, if the viscosity of the fracturing fluid continues to increase, although the leak-off situation continues to weaken, the friction resistance increases during hydraulic fracture expansion, making it difficult for the fracture to extend forward. Hence, 15 mPa s was considered to be the optimal fracturing fluid viscosity.

3.1.5 Perforation position

In on-site fracturing construction, a reasonable perforation location design has a direct impact on the fracture expansion and mining efficiency of CBM wells (Wang et al. 2019; Zhao et al. 2019). The perforation position refers to the distance between the perforation point and the upper horizontal interface of the coal seam. Figure 21 shows different perforation positions of 0.5, 1.25, 2.5 and 5.0 m, respectively. Numerical simulation of fracture extension under different perforation positions is performed. The simulated results are shown in Figs. 22 and 23.

Diagram of different perforation positions

Fracture morphologies under different perforation positions

Injection pressure curves under different perforation positions

Figure 22 shows that when the perforation location is 0.5 m away from the top of the coal seam, the formed fracture is similar to an H-shaped fracture. The upper and lower horizontal interface fractures of the H-shaped fracture are asymmetric. The vertical fracture height in the sandstone roof is not limited. When the distance between the perforation point and the upper horizontal interface of the coal seam increases from 1.25 to 2.5 m, the hydraulic fracture morphology changes from a cross-shaped fracture to a “士”-shaped fracture, and the hydraulic fracture in the coal seam is well developed. When the perforation location increases to 5 m away from the top of the coal seam, the hydraulic fracture morphology is dominated by vertical fractures, and the hydraulic fracture expands less in the coal seam. A large amount of fracturing fluid expands and filters out in the sandstone roof, triggering a serious decrease in the contribution to create hydraulic fractures.

Figure 23 shows the injection pressure curves under different perforation positions. As the perforation location increases from 0.5 to 5 m, the time required for the fracturing fluid to reach the interface increases, and the layer-penetrating pressure increases. When the hydraulic pressure is high enough, hydraulic fractures will penetrate the strata interface to propagate, thereby generating a vertical fracture. When the perforation location increases to 5 m away from the top of the coal seam, the hydraulic fracture morphology is dominated by vertical fractures, and the hydraulic fracture expands less in the coal seam.

Figure 24a shows that the fracture extension in the coal seam is relatively better when the distance between the perforation point and the upper horizontal interface is 0.5 m and 2.5 m. However, when the distance between the perforation position and the upper horizontal interface is 0.5 m, the vertical fracture length in the coal seam is 10.02 m, which is shorter than that of 2.5 m, and the horizontal fracture length in the upper horizontal fracture is 5.03 m, which is shorter than that of 2.5 m. Therefore, the optimal perforation position is selected at a distance between the perforation point and the upper horizontal interface of 2.5 m, which can obtain the maximum communication range of fractures in the coal seam.

Hydraulic fracture profiles under different perforation positions

3.2 Fracture morphology evolution law

According to the results of numerical simulations, it can be concluded that different parameter combinations may form different fracture morphologies, as illustrated in Fig. 25 and Table 3. When the perforation position is near the upper lithologic interface, it is prone to generate an H-shaped fracture. When the perforation location is far from the upper lithologic interface, it is likely to form a vertical fracture. When the coal seam permeability is lower (5–9 mD) and the minimum horizontal principal stress between layers and fracturing fluid injection rate are larger, it tends to form a “士”-shaped fracture. When the coal seam permeability is higher (9–30 mD) and the minimum horizontal principal stress between layers and fracturing fluid injection rate are smaller, it is apt to form a T-shaped fracture. When the coal seam permeability and minimum horizontal principal stress between layers and the perforation position are moderate, cross-shaped fractures will be easily generated. Therefore, this finding can provide a reference for hydraulic fracturing design.

Fracture morphology evolution under different parameters

4.1 Combination weight method

Several reliable methods have been widely applied to calculate the weights of factors affecting the shape parameters of hydraulic fractures. These methods have their own advantages and disadvantages. The combined weighting method combines several different methods to compensate for the shortcomings of each method and effectively export the advantages of each method. Among them, the weights obtained by the grey incidence analysis method take into account the impact of the changes in various factors on the parameters of hydraulic fractures. The AHP method establishes a relationship between various factors based on expert experience. The weight obtained by the entropy weight method is completely based on the relationship between the data.

4.1.1 Grey incidence analysis method

The basic idea of grey incidence analysis is to assess whether different data sequences are closely associated on the basis of the similarity of the geometric shapes of the reference series and several comparison series, reflecting the correlation degree of the curves (Liu et al. 2017). The calculation steps of grey incidence analysis are as follows.

4.1.1.1 Determination of the analysis sequences

Determine the reference sequence and the comparison sequence. A data sequence that reflects the characteristics of the behavior of the system is called a reference sequence. A comparison sequence is defined as a sequence that affects the system behavior.

The reference sequence (also known as the parent sequence) is

The comparison sequence (also called the child subsequence) is

4.1.1.2 Nondimensionalization of the variables

The data of each factor in the system may be different due to their different dimensions, and it is not convenient to compare or it is difficult to obtain a correct conclusion when comparing. Therefore, in the analysis of the gray incidence degree, nondimensional processing of the data is generally needed. There are two main processing methods: initial value and average value. The average value method was employed in this study.

4.1.1.3 Calculate the incidence coefficient

The so-called incidence coefficient is essentially the degree of difference in geometric shapes between curves. Therefore, the difference between the curves can be used as a measure of the incidence degree. The incidence coefficient \({\gamma }_{0i}\left(k\right)\) of each reference sequence and comparison sequence at each time (each point in the curve) (n = k) can be expressed as

where, \(\xi \in \left(\mathrm{0,1}\right)\) is the distinguishing coefficient; the smaller the \(\xi\) is, the greater the resolution, usually taking \(\xi =0.5\). \(\gamma \left({X}_{0},{X}_{i}\right)\) is the degree of grey incidence between \({X}_{0}\) and \({X}_{i}\), abbreviated as \({\gamma }_{0i}\). The k-point incidence coefficient \(\gamma \left({x}_{0}\left(k\right),{x}_{i}\left(k\right)\right)\) is abbreviated as \({\gamma }_{0i}\left(k\right)\).

4.1.1.4 Calculate the degree of grey incidence

Because the incidence coefficient is the correlation value between the comparison sequence and the reference sequence at each time (each point in the curve), the number of incidence coefficients is more than one, and the information is too scattered to facilitate overall comparison. Therefore, it is necessary to concentrate the incidence coefficients at each moment (each point in the curve) into one value, that is, to calculate the average value, as the quantitative expression of the incidence degree between the comparison series and the reference series, the calculation formula is as follows.

4.1.1.5 Sorting the incidence degree

The degree of grey incidence between factors is mainly depicted by the order of incidence degree, not just the size of incidence degree. Arrange the incidence degree of the m child subsequences to the same parent sequence in order of magnitude to form an association sequence, denoted as \(\left\{X\right\}\), which reflects the superiority relationship of each child subsequence to the parent sequence. If \({\gamma }_{01}>{\gamma }_{02}\), it represents that \(\left\{{X}_{1}\right\}\) is superior to \(\left\{{X}_{2}\right\}\) for the same parent sequence \(\left\{{X}_{0}\right\}\), that is, the reference sequence \(\left\{{X}_{0}\right\}\) is more similar to the comparison sequence \(\left\{{X}_{1}\right\}\).

4.1.1.6 Calculate the weight

After the incidence degree is normalized, the weight of the incidence degree of each comparison sequence is obtained.

where \({W}_{\mathrm{G}}(i)\) is the weight of the influencing factors calculated by the grey incidence analysis method.

4.1.2 Analytic hierarchy process (AHP) method

AHP is a decision analysis method that combines qualitative analysis with quantitative analysis methods to solve complex multiobjective problems (Process and Saaty 2013). This method uses the experience of experts to judge the relative importance of the standards that measure whether the objectives can be achieved. In this method, the importance of each factor needs to be subjectively judged first, and the discriminant result is represented by a numerical value. After processing the numerical value, the weight of each factor can be obtained. The calculation steps of the AHP method are as follows (Yi et al. 2020).

(1) Create a scale matrix. A is assumed to be the set of factors affecting hydraulic fracture parameters, including \({x}_{1},{x}_{2},\dots ,{x}_{m}\), a total of m factors, using the experience of experts and comparing the relative importance of elements to A, that is, the scale \({a}_{ij}\). \({a}_{ij}\) and \({a}_{ji}\) are reciprocals of each other, thus establishing a scale matrix \(\boldsymbol{A}=\left[{a}_{ij}\right]\). The corresponding relationship between scale assignment and element comparison importance is shown in Table 4. When \({a}_{ij}\) =2, 4, 6, 8, 1/2, 1/4, 1/6, 1/8, the corresponding importance level is between the levels of importance given in the table.

(2) Establishment of the attribute judgment matrix. The relative attribute is related to \({a}_{ij}\), so the attribute judgment matrix \(\boldsymbol{X}=\left[{a}_{ij}\right]\) can be established.

where, \({x}_{ij}\) is the relative attribute, \(j=\mathrm{1,2},\cdots ,m\); β is the correction coefficient of the judgment matrix, which is 1 in the text; and c is a constant.

(3) Calculate the weight of each factor.

where \({W}_{\mathrm{A}}(i)\) is the weight of each factor in AHP.

4.1.3 Entropy weight method

The entropy weight method is a method that determines the weight by the information provided by the observation values of various indicators. Entropy is a measure of the degree of disorder in the system. If the information entropy of an indicator is smaller, the amount of information provided by the indicator will be greater, its role in the evaluation will be greater, and the weight will be higher (Qiao et al. 2021). Therefore, when the entropy weight method calculates the weight, it only establishes the correlation between the parameters and does not consider the situation of the dependent variable, which is a purely theoretical weight calculation method. The calculation procedures of the entropy weight method are as follows (Xue et al. 2021).

1. (1)

   To avoid the difference in magnitude and unit between the data from affecting the readability of the data, the data are dimensionless.

2. (2)

   Calculate the proportion \({P}_{ij}\) of the index value of each simulation case under each influencing factor to all the index values under that factor.

   $${P}_{ik}=\frac{{b}_{ik}}{\sum_{k=1}^{n} {b}_{ik}}$$

   (10)

   where \({b}_{ik}\) is the dimensionless data of each factor and \({P}_{ik}\) is the proportion of the index value.

3. (3)

   Calculate the entropy value \({e}_{i}\) of each influencing factor.

   $${e}_{i}=-\frac{1}{\mathrm{ln}m}\sum_{k=1}^{n} {P}_{\mathrm{ik}}\times \mathrm{ln}\,{P}_{\mathrm{ik}}$$

   (11)

   where \({e}_{i}\) is the entropy value of influencing factor i, and the value is \(0\le {e}_{i}\le 1\).

4. (4)

   The larger the entropy value is, the smaller the difference of the factor evaluation value, and the smaller the effect it plays in the comprehensive evaluation. Therefore, the difference coefficient \({g}_{i}\) is introduced. The larger the difference coefficient is, the greater the influence of the factor, and the more important it is in the evaluation.

   $${g}_{i}=1-{e}_{i}$$

   (12)

   where \({g}_{i}\) is the difference coefficient.

5. (5)

   Calculate the weights based on the coefficient of variance.

   $${W}_{\mathrm{E}}(i)=\frac{{g}_{i}}{\sum_{i=1}^{m} {g}_{i}}$$

   (13)

   where \({W}_{\mathrm{E}}(i)\) is the weight of each factor in the entropy weight method.

Combining the influence indicators of fracture parameters, namely, multiplying the subjective and objective weights correspondingly, and normalizing the product, effectively avoids the influence of a single weight on the evaluation results.

where \(W(i)\) is the weight of the combination method.

4.2 Sensitivity analysis of fracture parameters

4.2.1 Grey incidence analysis of fracture shape parameters

The above experimental and numerical studies found that “干”-shaped fractures, T-shaped fractures, “士”-shaped fractures and cross fractures dominated by horizontal fractures were formed after roof fracturing, among which horizontal fractures were the main contributors to fracturing stimulation, Vertical seams are mainly responsible for the communication of coal seams. Moreover, the engineering scenarios of the Hancheng block also found that it is prone to form horizontal fractures in the lithological interfaces between the coal seam and the roof, and horizontal fractures play a major role in the fracturing effect (Zhao et al. 2019; Lyu et al. 2020).

Based on the numerical simulation cases in Sect. 3.1, the grey incidence analysis method was utilized to conduct a sensitivity analysis of the fracture shape parameters of roof fracturing. The fracture shape parameters of the upper horizontal fracture and vertical fracture are selected as the reference sequence, while their five influencing factors are selected as the comparison sequence. The calculated results are shown in Table 5. It is concluded that the incidence degrees of the five factors are all high, above 0.7, suggesting that the five factors are all closely related to the geometry of hydraulic fractures; additionally, different fracture parameters have different main controlling factors. For the half-length of the upper horizontal fracture, the fracturing fluid injection rate is the first important influencing factor, followed by the perforation location, and the fracturing fluid viscosity is the third. For the upper horizontal fracture width, the main controlling factors are the perforation position, fracturing fluid injection rate and coal seam permeability in order of importance. For the upper horizontal fracture height, the fracturing fluid injection rate, fracturing fluid viscosity and perforation position are the top three influencing factors. The perforation position has the strongest impact on the vertical fracture parameters. The minimum horizontal principal stress difference and fracturing fluid viscosity have a stronger influence on the vertical fracture length; fracturing fluid injection rate and viscosity are more important for the vertical fracture width; and fracturing fluid injection rate and minimum horizontal principal stress difference are more significant for the vertical fracture height. During fracture design, the desired fracture morphology can be obtained by adjusting the values of the corresponding main controlling factors.

4.2.2 Combination weight of fracture shape parameters

As listed in Table 6, the combined weight of engineering factors (perforation location, fracturing fluid injection rate and viscosity) is larger than the coal seam permeability and minimum horizontal principal stress difference between layers; perforation location has the strongest impact on the fracture shape parameters, followed by fluid injection rate and viscosity, while the influence of the coal seam permeability and minimum horizontal principal stress difference between layers is relatively weak. The weight ranking of the combination weight method is similar to those of the AHP method and grey incidence method, but different from that of the grey incidence method. The reliability of the weights of influencing factors of the fracture parameters obtained by the combined weight method is better than that of the single weight calculation method, which can effectively compensate for the solution deviation caused by the limited data samples or the difference in method principles of the single method.

The conclusions drawn from the physical and numerical simulation of the indirect fracturing in the roof of crushed soft coal seams are summarized as follows.

1. (1)

   “士”-shaped fractures, T-shaped fractures, cross fractures, H-shaped fractures, “干”-shaped fractures and vertical fractures are commonly generated during roof fracturing. The fracture morphologies can be transformed into each other under different geological and engineering conditions.

2. (2)

   Different parameter combinations can form different fracture morphologies. When the perforation position is near the upper lithologic interface, it is prone to generate H-shaped fractures. When the perforation location is far from the upper lithologic interface, it is likely to form vertical fractures. When the coal seam permeability is lower (5–9 mD) and the minimum horizontal principal stress between layers and fracturing fluid injection rate are larger, it tends to form “士”-shaped fractures. When the coal seam permeability is higher (9–30 mD) and the minimum horizontal principal stress between layers and fracturing fluid injection rate are smaller, it is apt to form T-shaped fractures. When the coal seam permeability and minimum horizontal principal stress between layers and the perforation position are moderate, cross-shaped fractures will be easily generated.

3. (3)

   For indirect fracturing in the roof of the No. 5 coal seam of the Hancheng block, the minimum horizontal principal stress difference between layers of 2 MPa, coal seam permeability of 15 mD, fracturing fluid injection rate of 3 m³/min, fracturing fluid viscosity of 15 mPa s, and perforation position of 2.5 m away from the upper interface of the coal seam were demonstrated to be the optimal parameters.

4. (4)

   The engineering factors of perforation location, fracturing fluid injection rate and viscosity have a strong impact on the fracture shape parameters. Different fracture shape parameters have different main controlling factors. For the half-length, width, and height of the upper horizontal fracture, the fracturing fluid injection rate, perforation position and fracturing fluid injection rate are the first important influencing factors, respectively. The perforation position, minimum horizontal principal stress difference, and fracturing fluid injection rate have the strongest impact on the length, width and height of the vertical fracture, respectively. During fracture design, the desired fracture morphology can be obtained by adjusting the values of the corresponding main controlling factors.

5. (5)

   In fracturing design, actual formation conditions should be fully considered, such as reservoir physical properties, in-situ stress distribution and sweet spot screening, and appropriate fracturing construction parameters, such as injection rate and viscosity of fracturing fluid and perforating parameters, should be optimized according to target fracture characteristics. In the on-site fracturing operations, real-time analysis of the pumping pressure curve, dynamic understanding of reservoir characteristics, and timely adjustment of operational parameters are required to ensure that the quantities of materials into the CBM well conform to the fracturing construction design, providing a reference for roof indirect fracturing design and construction of crushed soft coal seams.