Analysis of Crack Initiation and Propagation Thresholds of Inclined Cracks under High-Pressure Grouting in Ordovician Limestone

: Control of the grouting pressure within the critical grouting pressure for crack propagation in Ordovician limestone can not only ensure grout penetration length, but also prevent the risk of creating an artiﬁcial water channel. Based on the fracture mechanics theory, a formula was proposed to calculate the critical grouting pressure of mixed mode I-II cracks in Ordovician limestone. The necessary conditions for tilted crack opening, the rationality of the existing empirical value of the maximum allowable grouting pressure was investigated based on the mechanical model. The RFPA 2D -Flow numerical simulation software was used to evaluate the deduced theory. The research results show that the deduced theoretical calculation formula of the critical grouting pressure agrees with the numerical simulation results; when the mixed mode I-II fracture initiation occurs, the grouting pressure exceeds the perpendicular stress of the overlying rock; the greater the density of the overlying rock mass, the greater the value of grouting pressure for fracture initiation; when the side pressure coefﬁcient was ≥ 1, crack dip angle increased and the grouting pressure for fracture initiation tended to decrease; and the empirical grouting pressure at the maximum allowable grouting pressure is 2.0–2.5 p w , which will not cause propagation and failure of the existing crack.


Introduction
As the mining depth increases, coal mines in the North China coalfields are mainly exploiting Carboniferous Permian coal seams, lying above the Ordovician karst aquifer, which is a massive threat to the top ultra-close seam coal mining [1,2]. To prevent water inrush from the Ordovician karst aquifer, the grouting reinforcement technology of the upper part of Ordovician limestone is widely applied [3][4][5][6]. Due to the high water pressure in the Ordovician karst aquifer, the traditional grouting reinforcement of the coal floor cannot prevent the water inrush. Thus, several new drilling technologies have been developed to improve the grouting quality, including surface directional drilling, underground directional drilling, and radial jetting [5,7]. The slurry parameters in grouting holes, such as grouting pressure, grouting ratio, and grouting time, affect the grouting project's quality [8]. In recent years, many studies have been done on grouting materials' characteristics in cracks [9][10][11][12]. Li [13] studied the grout transport law from grouting holes in the Ordovician karst. Qiu [14] explored the characteristics of grouting reinforcement engineering practice in the Ordovician top. However, the above studies never focused on studying the critical grouting pressure in the Ordovician aquifer.
The main purpose of the grouting reinforcement of the Ordovician top is to achieve sealing and reinforcement of the karst fracture network, thus acting as an anti-leakage plug [8,15]. Normally, as the grouting pressure increases, the grout penetration length is increased, which helps improve the strength of the skeleton and make the rock and soil mass more compact and solidified. When the grouting pressure exceeding the fracture initiation, grouting pressure may lead to hydraulic fracture propagation [16], resulting in cleavage and the formation of artificial water-conducting crack channels, causing secondary damage. The best quality of grouting projects can be achieved with a suitable penetration length by controlling the grouting pressure in grouting projects. Few reliable models can be applied to calculate the critical grouting pressure for inclined fissures in the Ordovician limestone [17], resulting in the practical value of the maximum allowable grouting pressure 2.0-2.5 p w (water pressure) become the unique, reliable value in grouting engineering in North China [8]. Until now, the reliability of this empirical value has not been demonstrated.
The Ordovician limestone open cracks caused by grout pressure damage by belong to the fluid high-pressure fracturing and propagation type. Due to the complex loading conditions and lithologic heterogeneity of the real rock, simple mode I (opening), II (sliding), and III (tearing) cracks rarely exist [18], instead mixed-mode I-II cracks are mainly found [19]. Several mixed-mode crack criteria were developed to calculate fracture initiation's critical hydraulic pressure [18,[20][21][22]. Erdogan and Shi [23] analyzed for the first time the mixedmode I-II crack propagation. Since then, the minimum strain energy density criterion [24], the maximum potential energy release-rate criterion [25], the maximum tangential strain criterion [26], and the T-criterion [27] have been proposed, respectively. The above criteria often involve sequential procedures, bringing difficulties for the critical hydrostatic pressure calculation when applied to crack initiation analysis. Zhou [28] created an empirical criterion for crack initiation and propagation through rock mechanics tests. The empirical standard is supported by other observations [29,30], and it is widely used because of its simple but useful expressions. However, none of the above models has been introduced to solve the limestone crack initiation problem.
In general, the scale of cracks in limestone ranges from millimeters to kilometers. It is difficult to simulate the grouting process in meter-scale cracks with different crack lengths, crack angles and crack stress state in laboratory experiments [16]. To better verify the created fracture criterion, numerical methods are applied in the research of hydraulic industries. Finite Difference Method (FDE), Finite Discrete Element Method (FDEM), Discrete Element Method (DEM) and Half Analytical Element Method (HAEM) have been used to model the composite crack initiation and propagation during hydraulic fracturing [31,32]. The above numerical methods are useful, but there are still some shortcomings in solving the problems of nonlinearity and discontinuities of rocks. To solve the above problems, the Realistic Failure Process Analysis (RFPA) method was developed by Tang [33].
The purposes of this study were: (1) to derive the crack initiation grouting pressure formula based on Zhou's empirical criterion for calculating the critical grouting pressure of mixed mode I-II cracks in Ordovician limestone; (2) to examine the destruction mode of crack propagation under the high-pressure grouting conditions; (3) to analyse the influence of factors such as crack dip angle, overlying rock density, and lateral pressure coefficient on the fracture initiation and (4) to verify the rationality of the derived theory and the existing empirical value of the maximum allowable grouting pressure (2.0-2.5 p w ) by using the RFPA2D-Flow simulation software.

Determining the Crack Initiation Grouting Pressure Using Zhou's Criterion
In natural saturated cracks, the hydrostatic pressure in the crack is in equilibrium with the force exerted by the surrounding rock; a rock with a crack can still be considered as an intact rock, which can be analysed by rock mechanics, and a single inclined crack in a rock can be considered as a stress element [34]. The stress analysis of a single inclined crack under the grouting is shown in Figure 1, with σ 1 as the vertical principal stress, σ 3 as the horizontal principal stress, and α as the crack dip angle. The angle between the inclined crack and the maximum principal stress is represented by α; the crack dip angle, by β; the positive stress in the normal direction of the crack plane, by σ n ; shear stress, by τ α ; and pore water pressure, by p.

Determining the Crack Initiation Grouting Pressure Using Zhou's Criterion
In natural saturated cracks, the hydrostatic pressure in the crack is in equilibrium with the force exerted by the surrounding rock; a rock with a crack can still be considered as an intact rock, which can be analysed by rock mechanics, and a single inclined crack in a rock can be considered as a stress element [34]. The stress analysis of a single inclined crack under the grouting is shown in Figure 1, with σ1 as the vertical principal stress, σ3 as the horizontal principal stress, and α as the crack dip angle. The angle between the inclined crack and the maximum principal stress is represented by α; the crack dip angle, by β; the positive stress in the normal direction of the crack plane, by σn; shear stress, by τα; and pore water pressure, by p. By using Roylance's [35] coordinate conversion formula, the positive and shear stresses at the cracked surface after coordinate conversion conform to the following relational formula: Fracture mechanics specifies that tension is positive, and pressure is negative, which is opposite to rock mechanics. Hence, a negative sign is added in front of Equations (1) and (2). The mixed mode I-II crack propagation can be divided into tensile-shear and pressure-shear modes [36]. When the normal positive stress at the crack surface is tensile, the mixed mode I-II crack can be classified as tensile-shear mode, and when the normal positive stress at the crack surface is compressive, it belongs to pressure-shear mode.
Grout slurry can only transport in Ordovician limestone open cracks [15], of which the faces are not in contact, and the extension loading tends to be tensile-shear mode [37]. When the cracks slip, there is an obvious shear dilation at the crack tip, but the stress on the cracks will stop this expansion. Due to the complexity of the compressive-shear fracture mechanism of rocks, there is no well-established and generally accepted criteria for mixed-mode I-II cracks failure. Zhou's empirical failure criterion method is one of the more practical ways of solving crack in solving the problem in fracture mechanics. Since By using Roylance's [35] coordinate conversion formula, the positive and shear stresses at the cracked surface after coordinate conversion conform to the following relational formula: Fracture mechanics specifies that tension is positive, and pressure is negative, which is opposite to rock mechanics. Hence, a negative sign is added in front of Equations (1) and (2). The mixed mode I-II crack propagation can be divided into tensile-shear and pressure-shear modes [36]. When the normal positive stress at the crack surface is tensile, the mixed mode I-II crack can be classified as tensile-shear mode, and when the normal positive stress at the crack surface is compressive, it belongs to pressure-shear mode.
Grout slurry can only transport in Ordovician limestone open cracks [15], of which the faces are not in contact, and the extension loading tends to be tensile-shear mode [37]. When the cracks slip, there is an obvious shear dilation at the crack tip, but the stress on the cracks will stop this expansion. Due to the complexity of the compressive-shear fracture mechanism of rocks, there is no well-established and generally accepted criteria for mixed-mode I-II cracks failure. Zhou's empirical failure criterion method is one of the more practical ways of solving crack in solving the problem in fracture mechanics. Since other failure criteria assume that the rock is homogeneous and anisotropy, but the Zhou's empirical failure criterion could ignore the above hypothesis. For this reason, the empirical failure criterion was selected to describe the crack initiation in the in Ordovician limestone. The Zhou's empirical failure criterion calculated the crack initiation by the following formula [28,38]: where K IC is mode I fracture toughness, MPa·m 1/2 ; K I is crack mode I stress intensity factor, MPa·m 1/2 ; and K II is crack mode II stress intensity factor, MPa·m 1/2 . Based on fracture mechanics, the mode I and II stress intensity factors are calculated as follows [39]: where a is the half-length of the crack m. Substituting Equations (1) and (2) into (4) and (5), respectively, yields the following formula for calculating the critical grouting pressure for mixed-mode I-II crack initiation: By defining the side pressure coefficient as k = σ 3 σ 1 , the Equation (6) becomes: The critical grouting pressure for mixed-mode I-II crack is obviously related to overburden stress, fracture dip, fracture size and the side pressure coefficient. The effects of overburden stress, fracture dip, and lateral pressure coefficient on the crack initiation grouting pressure are discussed later in Section 4.

Mechanical Criteria for Crack Opening
During the grouting process of the limestone cracks, when the grouting pressure is low, the grout flow rate is small, and infiltration and filling dominate. When the grouting pressure is high, the fracture surface will be deformed to a certain extent, also known as crack deformation. Crack deformation is mainly an elastic deformation.
Some scholars [40] argued that grout shows a radiation pattern when it diffuses in sheet cracks, the grouting pressure on the crack surface is not equal, and the grout pressure gradually decreases with the increase in diffusion radius; that is, P 0 ≥ p, as shown in Figure 2a. The width of the crack propagation decreases with the increase in the diffusion radius. As the crack rock mass mainly undergoes elastic deformation, the grout-rock interface deformation process can be simplified to a semi-infinite space elastomer pressure model, as shown in Figure 2b. An oxy-rectangular plane coordinate system is established at the grout-rock interface (Figure 2b), and only the additional stress generated by the pressure of the grout in the crack after overcoming the stresses on the surrounding rock can cause deformation of the grout-rock interface. The additional stress formula should meet the following condition: An oxy-rectangular plane coordinate system is established at the grout-rock interface (Figure 2b), and only the additional stress generated by the pressure of the grout in the crack after overcoming the stresses on the surrounding rock can cause deformation of the grout-rock interface. The additional stress formula should meet the following condition: where σ α is the additional stress perpendicular to the z direction of the oxy plane, Mpa, and p is the grout pressure in the crack, Mpa. The simplified Equation (8) shows that the grout pressure in the crack must meet the following condition if the crack is to be opened or deformed: Substituting the Equation (2) into (9), the necessary condition for the deformation of the fracture of the rock mass is In near-surface or shallow strata, the horizontal stress is generally greater than the vertical stress [41]; that is, σ 3 ≥ σ 1 . Then, Equation (10) is transformed to from Equation (11), σ 3 ≥ σ 1 , where σ 1 is the self-weight pressure of the surrounding rock, which can be calculated from the following equation: where γ R is the average bulk density of the rock overlying the grouting section (N/m 3 ) and h is the thickness of the stratum above the grouting section (m). Substituting Equations (10) and (12) into (11), the necessary condition for the deformation of the fracture surface during the grouting of the fractured rock mass is

Analysis of the Rationality of the Maximum Allowable Grouting Pressure
The maximum allowable grouting pressure refers to the maximum grouting pressure allowed by the boundary conditions during the grouting process. The control of grouting pressure is the key success factor to the grouting project. In the Ordovician limestone grouting, the grout must be controlled within the maximum allowable grouting pressure range. The geological conditions, method of grouting, and concentration of the grout all affect the choice of the grouting pressure, and a more reasonable method is to determine the grouting pressure through field tests. Generally, the following empirical formula is used for calculation [8]: where p max is the maximum allowable grouting pressure at the crack inlet (MPa) and p w is the groundwater hydrostatic pressure at the grouting segment (MPa). The groundwater static water pressure in the grouting segment is calculated as: where γ w is the water capacity and h w is the hydrostatic water level (m). When h w = h, the hydrostatic pressure p w reaches its maximum value, the maximum allowable grouting pressure is taken according to the empirical value, and the maximum value is taken as p max = 2.5 p w . Thus: As the average weight capacity of the rock mass γ R = (2.5~3)γ w , we can obtain σ 3 = (2.5~3)γ w h, which shows that the maximum grouting pressure is difficult to exceed the self-weight stress of the overlying rock mass. According to Equation (16), the necessary condition for deformation is a grouting pressure greater than the self-weight stress of the overlying rock. Thus, in deep grouting projects, as long as the grouting pressure does not exceed the maximum allowable grouting pressure, the crack grout-rock interface will not deform; that is, the elastic deformation of the crack rock mass is negligible in the crack grouting project, and the crack width b can be considered as a constant in the grouting project. Whether p max = 2.5 p w is appropriate for propagation when the overlying rock layer has a bulk density of <2500 kg/m 3 , is verified later by numerical simulation.

Study Area
The Xingtai coalfield is located in Xingtai City of Hebei Province, China ( Figure 3). It is a typical North China coalfield, where the coal seams are in the Permo-Carboniferous. There are 10 mines in the Xingtai coalfield, and the mining depth is more than 600 m. The stratigraphy of the Xingtai coalfield mainly includes the Ordovician system (O), Carboniferous system (C), Permian system (P), Paleogene system (E), and Cenozoic Quaternary (Q). The No. 9 coal seam is one of three main mineable coal seams, were threatened by the main aquifers in the Ordovician limestone with 7.0 Mpa water pressure. The distance from the No. 9 coal seam floor to the top of the Ordovician strata is about 40-50 m [4,5,42].

Numerical Model
To better evaluate the deduced theory for crack initiation in limestone and verify the rationality of the existing empirical value of the maximum allowable grouting pressure. The RFPA 2D -Flow numerical simulation software was used for its special features that consider the nonlinearity, non-homogeneity and the anisotropic properties of rocks [33]. Besides, several reports have shown that limestone has the characteristics elastic-brittle,

Numerical Model
To better evaluate the deduced theory for crack initiation in limestone and verify the rationality of the existing empirical value of the maximum allowable grouting pressure. The RFPA 2D -Flow numerical simulation software was used for its special features that consider the nonlinearity, non-homogeneity and the anisotropic properties of rocks [33]. Besides, several reports have shown that limestone has the characteristics elastic-brittle, or brittle-plastic [16,43]. The RFPA method chooses Modified Mohr-Coulomb as the damage criterion [19,44].
A computational model is established, as shown in Figure 4. The model size is 0.6 × 1 m, the number of grid divisions is 150 × 250 = 375,000, the model crack length is 2a = 0.2 m, and the crack width is 2 mm, suggesting that the crack is narrow and suitable for grouting [45].  The numerical simulation is divided into two groups: limestone crack with different dip angles (0°, 15°, 30°, 45°, 60°, and 75°) under unchanged external pressure conditions and cracks with different overburden pressures by changing the average density of the overlying rock. A total of 10 scenarios were simulated in this numerical simulation experiment, and the details are shown in Table 1.

Model Setup
(1) Material parameters Grouting projects are mostly for the limestone area management in North China [5,7], and the crack medium for grouting transformation and reinforcement is mostly limestone mass; therefore, the numerical simulation adopts the relevant parameters of limestone [46][47][48]. The mechanical and permeability parameters are shown in Table 2. Owing to the inhomogeneity of the rock mass, the inhomogeneity degree m = 5 is taken; the simulated rock mass is 0 flow boundary, and the grout density is 1.216 kg/m 3 .  Table 1.

Model Setup
(1) Material parameters Grouting projects are mostly for the limestone area management in North China [5,7], and the crack medium for grouting transformation and reinforcement is mostly limestone mass; therefore, the numerical simulation adopts the relevant parameters of limestone [46][47][48]. The mechanical and permeability parameters are shown in Table 2.
Owing to the inhomogeneity of the rock mass, the inhomogeneity degree m = 5 is taken; the simulated rock mass is 0 flow boundary, and the grout density is 1.216 kg/m 3 . (2) Boundary conditions The average mining depth of coal mines in Xingtai coalfield is approximately 0.6 km [17,49], the simulation set buried depth of the limestone cracks to be 600 m, and average density of the overburden rock to be 2300, 2500, 2700, and 2900 kg/m 3 , corresponding to the overburden pressures at 13.8, 15, 16.2, and 17.4 Mpa, respectively. The horizontal stress is set according to the relationship of burial depth. Previous studies showed that the relationship between horizontal stress and depth is a linear function, and maximum horizontal principal stress and minimum horizontal stress have the following relationships with burial depth [41,50,51]: where σ H is the maximum horizontal principal stress and σ h is the minimum horizontal stress. Taking the burial depth H = 600 m into the formula, the average horizontal stress is calculated to be 16.2 Mpa.

(3) Solution Settings
During the simulation process, gravity and initial stress on crack propagation in the grouting process are considered. Before the grouting pressure is applied, 6 steps have to be run to achieve redistribution of the initial stress field. The increment of grouting pressure ∆p is 0.1 Mpa (at a grout column height of 8.224 m). The total simulation steps are determined according to the specific simulation conditions, and the fluid-structure coupling is used to solve the problem.

Effect of Overburden Pressures on Fracture Initiation Grouting Pressure
In the RFPA 2D -Flow system, the location of fracture initiation and number of fracture initiation calculation steps can be found from the acoustic emission map. The number of fracture initiation calculation steps for different overlying rock densities can be obtained by exporting the acoustic emission data from the numerical simulation system and plotting the acoustic emission energy accumulation map. Figure 5 shows the number of crack initiation steps for rock failure under four different overlying rock densities with a dip angle of 45 • . It can be seen from the figure that the overburden pressures, corresponding to different overburden density, has a certain influence on the propagation and fracture initiation of the crack, which indicates that when the crack is located at a certain crack depth, the higher the average density of the overburden rock, the more calculation steps According to the fracture initiation steps obtained from the acoustic emission data, the fracture initiation grouting pressures under different average overlying rock densities are calculated and compared with the calculated values of the fracture initiation grouting pressure theory (Equation (6)). The data in Table 3 shows that (1) as the average density of the overburden increases, representing the increase of overburden stress, the fracture initiation grouting pressure value of the crack instability increases, and the fracture initiation grouting pressure is greater than the self-weight stress of the overlying rock, (2) under the same overburden density, the overall error between the calculated and numerical simulation values is <6%, and the maximum error value is 1.263 Mpa. This indicates that the calculated value of the fracture initiation grouting pressure is close to the numerical simulation value, and the limestone calculation formula has certain reliability when applied to calculate the grouting pressure value.  Figure 6 shows the comparison curves of the theoretical and numerical simulation values of the fracture initiation grouting pressure with the empirical value of the maximum allowable grouting pressure under different average density, representing different overburden pressure. It can be seen from Figure 6 that the grouting initiation pressure of the limestone crack is generally >3 times the hydrostatic pressure (3 pw). When the average density is 2300 kg/m 3 (overburden stress is 1.8 Mpa), the theoretical calculation value is 3.35 times the hydrostatic pressure (3.35 pw), and the numerical simulation value is 3.38 times the hydrostatic pressure (3.38 pw). As the overburden stress (average density) of the overlying rock increases, the theoretical and numerical simulation values of the grouting pressure both increase. The above results show that with the overburden density of 2300-3400 kgm −3 , when the grouting pressure is 2.5 pw, the crack propagation will not occur at According to the fracture initiation steps obtained from the acoustic emission data, the fracture initiation grouting pressures under different average overlying rock densities are calculated and compared with the calculated values of the fracture initiation grouting pressure theory (Equation (6)). The data in Table 3 shows that (1) as the average density of the overburden increases, representing the increase of overburden stress, the fracture initiation grouting pressure value of the crack instability increases, and the fracture initiation grouting pressure is greater than the self-weight stress of the overlying rock, (2) under the same overburden density, the overall error between the calculated and numerical simulation values is <6%, and the maximum error value is 1.263 Mpa. This indicates that the calculated value of the fracture initiation grouting pressure is close to the numerical simulation value, and the limestone calculation formula has certain reliability when applied to calculate the grouting pressure value.  Figure 6 shows the comparison curves of the theoretical and numerical simulation values of the fracture initiation grouting pressure with the empirical value of the maximum allowable grouting pressure under different average density, representing different overburden pressure. It can be seen from Figure 6 that the grouting initiation pressure of the limestone crack is generally >3 times the hydrostatic pressure (3 p w ). When the average density is 2300 kg/m 3 (overburden stress is 1.8 Mpa), the theoretical calculation value is 3.35 times the hydrostatic pressure (3.35 p w ), and the numerical simulation value is 3.38 times the hydrostatic pressure (3.38 p w ). As the overburden stress (average density) of the overlying rock increases, the theoretical and numerical simulation values of the grouting pressure both increase. The above results show that with the overburden density of 2300-3400 kgm −3 , when the grouting pressure is 2.5 p w , the crack propagation will not occur at the tip of natural limestone crack, which means the maximum allowable grouting pressure, 2.0-2.5 p w , is reasonable. the tip of natural limestone crack, which means the maximum allowable grouting pressure, 2.0-2.5 pw, is reasonable.  Figure 7 shows the number of calculation steps of fracture initiation in the numerical simulation acoustic emission data for cracks with dip angles of 0°, 15°, 30°, 45°, 60°, and 75°. As shown in Figure 6, under fixed surrounding rock conditions, the number of steps required for crack initiation tends to decrease as the crack dip angle increases, from 166 steps at 0° to 149 steps at 75°.      Figure 7 shows the number of calculation steps of fracture initiation in the numerical simulation acoustic emission data for cracks with dip angles of 0°, 15°, 30°, 45°, 60°, and 75°. As shown in Figure 6, under fixed surrounding rock conditions, the number of steps required for crack initiation tends to decrease as the crack dip angle increases, from 166 steps at 0° to 149 steps at 75°.   Table 4 shows the comparison of the values of the crack initiation grouting pressure with the calculated value (Equation (6)) and the numerical simulation results under different dip angle conditions with the side pressure coefficient k = 1. From Table 3, we can see that the overall error between the simulated and calculated values is <15%, and the maximum relative error is 13.82%, which shows that under different crack angle conditions, the calculation of the fracture initiation grouting pressure for mixed mode I-II crack is reliable.  Figure 8 shows the comparison curve between the calculated and numerical simulation values of crack initiation grouting pressure and the empirical value of the maximum allowable grouting pressure for different dip angle conditions.  Table 4 shows the comparison of the values of the crack initiation grouting pressure with the calculated value (Equation (6)) and the numerical simulation results under different dip angle conditions with the side pressure coefficient k = 1. From Table 3, we can see that the overall error between the simulated and calculated values is <15%, and the maximum relative error is 13.82%, which shows that under different crack angle conditions, the calculation of the fracture initiation grouting pressure for mixed mode I-II crack is reliable.  Figure 8 shows the comparison curve between the calculated and numerical simulation values of crack initiation grouting pressure and the empirical value of the maximum allowable grouting pressure for different dip angle conditions. From Figure 8, we can see that the fracture initiation grouting pressures are all greater than the weight stress of the overlying rock, and the propagation fracture initiation pressures of the crack are generally greater than three times the hydrostatic pressure (3 pw), which indicates that when the grouting pressure is 2.5 pw, the natural limestone cracks will not initiate and expand. As the crack dip angle increases, the numerically simulated value of the fracture initiation grouting pressure shows a decreasing trend, while the theoretically calculated value is a constant value 3.35 times the hydrostatic pressure (3.35 pw). The possible reason for this phenomenon is that the type I or II tensile-shear composite crack instability theory assumes that the crack is in a homogeneous rock mass; that is, the side pressure coefficient k = 1. At this time, the influence of each direction of the surrounding rock on the crack is equal, the shear stress at the crack surface is τα = 0 Mpa, and the value of the fracture initiation grouting pressure is not related to the dip angle of the crack. The RFPA 2D -Flow numerical simulation software takes into account the inhomogeneity of  From Figure 8, we can see that the fracture initiation grouting pressures are all greater than the weight stress of the overlying rock, and the propagation fracture initiation pressures of the crack are generally greater than three times the hydrostatic pressure (3 p w ), which indicates that when the grouting pressure is 2.5 p w , the natural limestone cracks will not initiate and expand. As the crack dip angle increases, the numerically simulated value of the fracture initiation grouting pressure shows a decreasing trend, while the theoretically calculated value is a constant value 3.35 times the hydrostatic pressure (3.35 p w ). The possible reason for this phenomenon is that the type I or II tensile-shear composite crack instability theory assumes that the crack is in a homogeneous rock mass; that is, the side pressure coefficient k = 1. At this time, the influence of each direction of the surrounding rock on the crack is equal, the shear stress at the crack surface is τ α = 0 Mpa, and the value of the fracture initiation grouting pressure is not related to the dip angle of the crack. The RFPA 2D -Flow numerical simulation software takes into account the inhomogeneity of lithology, and obey a certain distribution, e.g., Normal Distribution, Weibull Distribution or Equal Distribution. In our numerical simulation, the density uniformity of rock mass is set at 5.

Influence of the Fracture Dip Angle on the Fracture Initiation Grouting Pressure
The lateral pressure coefficient = 1 of the cracked surrounding rock directly leads to the correlation of the fracture initiation grouting pressure with the crack dip angle. The results of the studies of Tang et al. [18] and Li et al. [22] showed that when the lateral pressure factor was ≥1, the fracture initiation grouting pressure tended to decrease with the increase in crack dip angle. The results were consistent with the numerical simulation results, indicating that the RFPA 2D -Flow numerical simulation results are reliable.

Conclusions
Based on the Zhou's empirical rock failure criterion, a formula for calculating the critical grouting pressure for the mixed mode I-II crack was deduced, taking into account the influence of factors such as side pressure coefficient, fracture dip angle, and overburden stress. The theoretical formula is verified by numerical simulation. The calculation results were relatively consistent with the numerical simulation results, indicating that the deduced theoretical formula has a certain reference value.
The overburden stress on the crack and the dip angle of the crack influence the real limestone crack initiation grouting pressure; that is, the greater the overburden stress, the greater the value of fracture initiation grouting pressure. When the lateral pressure coefficient was ≥1, with the increase of the dip angle of the crack, the value of fracture initiation grouting pressure tended to decrease.
The numerical simulation and theoretical values of crack initiation showed that the fracture initiation pressure for crack propagation was greater than the vertical stress of the overlying surrounding limestone. The crack grouting initiation pressure was 2.5 times higher than the hydrostatic pressure, which is approximately 3.5-4 times higher than the hydrostatic pressure, indicating that fracture initiation will not occur when the maximum grouting pressure is 2.0-2.5 p w . The empirical value of maximum limestone grouting is reasonable.
The first limitation of this study is that the deduced formula is based on Zhou's empirical rock failure criteria, lacking mechanical analysis and comparison with other rock failure criterion, such as the minimum strain energy density criterion, the maximum potential energy release-rate criterion, the maximum tangential strain criterion, and the T-criterion. Moreover, which rock failure criterion can provide better the critical grouting pressure values need to be debated in future research. The second limitation of this study is that the deduced formula is only verified by the RFPA2D-Flow numerical simulation software, lacking experimental data support. Since there is always lacking indoor largescale grouting experiments, it is necessary to evaluate the deduced formula for the grouting initiation pressure using other approaches. One of the approaches is doing grouting tests in limestone cracks in the field using controlled trials.