Model tests and numerical simulations on hydraulic fracturing and failure mechanism of rock landslides

Hydraulic fracturing (HF) is an important reason for inducing rock landslides. The impacts of the angle between the transverse pressure stress and the initial fracture, the transverse pressure stress, the length and aperture of the initial fracture on the critical water pressure stress of HF were analyzed based on the model tests of cement mortar specimens. The HF numerical simulations were carried out and the fracture aperture, pore pressure, strain and stiffness degradation index were analyzed, and the HF mechanism of rock was revealed. Taking a rock slope along G205 in China as an example, the landslide stability under high groundwater pressure was studied by numerical simulations, and the whole processes of the expansion, propagation, penetration and the rock slope failure were reproduced. The results showed that when HF occurred on the cement mortar specimens, a large amount of water seeped from the fracture and made a dull sound of failure, the water pressure on the fracture plane decreased rapidly, but not completely penetrated, and there was still a certain residual strength of the specimen. The HF of rock is a quasi-brittle failure, including static stage, micro fracture propagation stage and macro fracture formation stage. There are 3 dangerous rocks on the rock slope, among which WYT3 is in a stable state under the natural condition and in an unstable state under the rainstorm condition and seismic condition. The fracture propagation of WYT3 under the action of high groundwater pressure experienced slow development stage, rapid development stage and penetration stage, among which the slow development stage lasted the longest, and the slope failure occurred immediately after the fracture penetrated.


Introduction
Rock landslides are serious natural hazards that have the potential to cause huge human casualties and economic losses (Yu et al. 2019). Globally, rock landslides cause approximately 300 deaths per year and a property damage of about 1 billion dollars (Zhan et al. 2019;Liao et al. 2019). According to Raja et al. (2017) and Zhang et al. (2020a, b, c, d), rock landslides account for 3% of all deaths caused by natural disasters. Hydraulic fracturing (HF) is a physical phenomenon of the initiation, propagation and penetration of micro fracture in rock under high water pressure until generate macro fracture and causes slope failure (Zhan and Cen 2007;Wang et al. 2020;Lyu et al. 2020). HF treatment is an efficient stimulation method in the development of unconventional hydrocarbon reservoirs with low permeability such as gas and oil shales, tight gas sands, coal bed methane, etc. (Dehghan 2020), however, HF is also an important reason for inducing structural failure such as tunnels, dams and rock landslides. For example, the Teton dam gushing in the United States, the Hyttejuvet dam leakage in Norway, the Jinping power plant gushing in China and the Valion landslide in Italy are catastrophic events caused by HF (Baghbanan and Jing 2007;Liu et al. 2018).
The essential reason for HF of rock is that the stress intensity factors at the fracture tip exceed the fracture toughness (Park et al. 2002). The studies of HF mainly adopt the following methods: in-situ test, model test, numerical simulation and analytical analysis. Among them, model tests can be used to study the splitting mode and the relationship between the critical water pressure stress and the confining pressure of rock specimens, thick-walled cylinder specimens and cement mortar specimens, which is the most intuitive method to reveal the HF mechanism of rock (Xu et al. 2017;Hou et al. 2019;Sun et al. 2020a, b;Zhuang et al. 2022), for example Liu et al. (2018) proved that the critical water pressure stress of HF was positively correlated with the material strength and the axial pressure stress, and negatively correlated with the length and the aperture of the initial fracture through model tests. Hou et al. (2019) simulated HF in different lithological rocks in the horizontal drilling by using the true physical model experiment and large rock specimens and carried out the real-time dynamic monitoring by adding tracer and then did post-fracturing cutting. Zhang et al. (2019a, b) analyzed the effects of different injection orifice positions on HF of shales, the results showed that the bedding plane was the hydraulic path of liquid flow, and multi-directional injection orifices could cause the development of complex fractures. Liu et al. (2019) performed a series of HF physical simulation experiments under true triaxial stress states, and acoustic emission devices were used to monitor real-time fracture propagation. Moreover, the key influencing factors and the typical fracturing curve characteristics of vug-HF interaction were studied and analyzed. Diaz et al. (2020) investigated the evolution of Acoustic Emission (AE) activity, hydraulic energy, fracture area and fracture development for cyclic injection schemes with different cycle Time Durations (TD) at a repeated maximum cyclic pressure through laboratory experiments. Hu et al. (2020) used model tests to analyze the HF deformation characteristics of dams under different water pressures and initial fracture lengths, the results showed that with the increases in water pressure and initial fracture length, the dam deformation increased, and the turn-to-different feature point of the whole displacement curve in advanced. Due to the disadvantages of HF model tests, such as the difficulties in sealing water and the difficulties in accurately measuring the pressure in the fracture and the strain at the fracture tip, the applications of HF model tests are limited to a certain extent. It is often necessary to verify the model test results with the corresponding numerical simulation results (Wu et al. 2017;Zou 2021).
The numerical model of HF simulation includes equivalent continuum medium model and discrete fracture medium model. Among them, the equivalent continuum medium model represents the fracture propagation by damages and establishes the relationship between the damage variables and permeability coefficients, while the discrete fracture medium model assumes that water exists and moves in rocks as the form of "fracture flow" (Hadjigeorgiou et al. 2009;Shen et al. 2019;Liu et al. 2021;Zhang et al. 2022). Wu et al. (2014) carried out the fluid-solid coupling simulations of the slope rainfall erosion process, and obtained important parameters such as porosity, particle trajectory and flow velocity when simulating large particle deformation, and the slope erosion degree and the distribution of erosion capacity during rainfall were revealed. Sun et al. (2022) established an innovative hydraulic aperture prediction equation for 2D fractures based on the mechanical aperture and fracture morphology, and 6 parameters were firstly selected to predict the hydraulic aperture, including 3 mechanical aperture parameters and 3 fracture morphology parameters. Yu et al. (2020) established a numerical model to study the propagation of non-planar fractures in stratified rock by using extended Finite Element Method (FEM). Xie et al. (2020) developed a fractal fluid flow model based on the anisotropic fractal characteristics of natural fracture planes and built 3D printed fracture models by using the morphological information of an actual tensile and shear fracture plane of a rock. Liu et al. (2020a, b) analyzed the dynamic process of the Xinmo landslide by using the PFC software, and the micro parameters of the rock landslide were calibrated with uniaxial compression and landslide deposition simulations. Chen et al. (2020) studied the fracture behavior of single-and multi-stage HF under varying excavation stress conditions by using a flow-coupled Discrete Element Method (DEM), and conventional HF theories and analytical solutions for excavation stress redistribution were used to verify the numerical model. The HF numerical simulations based on DEM are prone to the shortcomings of non-convergence at model boundary, difficult to conduct macro-scale simulation and long-time consuming, while FEM has advantages in solving a large range of flow field, so FEM is adopted to carry out HF numerical simulations (Xu et al. 2013;Ni et al. 2015;Alneasan et al. 2019;Zhang et al. 2020a, b, c, d).
The rock landslide failure requires not only mechanical conditions but also weak fracture planes that satisfy certain hydraulic conditions. The formation of weak fracture planes is closely related to the development process of slope fractures. When groundwater seeps into the fractures, the water pressure increases and the fractures continue to expand until they penetrate, resulting in an increase in the permeability of rock slopes and a faster seepage rate, eventually lead to the rock slope failure (Zhou et al. 2017;Zhang et al. 2019a, b;Qiu et al. 2020;Zhang et al. 2020a, b, c, d). Aiming at the complex coupling relationship between the fracture propagation and water seepage of rock, the influencing factors of HF were analyzed based on model tests and numerical simulations, and the HF mechanism was revealed. Then, numerical simulations were used to reproduce the whole processes of fracture propagation, penetration until failure of a rock slope along Letuan-Qingshiguan section of G205 in China. The research contents include HF model tests of cement mortar specimens; HF numerical simulations of cement mortar specimens; HF mechanism analysis of rock; numerical simulations of rock slope failure under high groundwater pressure.

Model preparation
Natural rock is rich in joints, so it is difficult to prepare a large number of specimens with similar size and physical properties, and artificial fracture is easy to cause damages. Thickwalled cylinder specimens are not difficult to make artificial fractures and hydraulic seal, but the relationship between the critical water pressure stress of HF and the length and aperture of the initial fracture cannot be measured. Solidified cement mortar has similar physical properties to rock and is easy to make artificial fracture. In this research, cement mortar specimens were selected for HF model tests (Zhan and Cen 2007;Liu et al. 2018). In order to facilitate the comparison between research results, scholars all over the world select the size of the specimen as 500 mm × 200 mm × 200 mm when conducting similar studies (Hadjigeorgiou et al. 2009;Liu et al. 2019;Shi et al. 2020). The mass ratio of the specimen was: cement: sand: water = 1.00: 4.93: 1.16. When pouring the specimen, a steel sheet of 100 mm width was embedded in the center of the mold to form the initial fracture. After 28 days of curing, the surface defects of the specimen were filled and polished with sandpaper. After testing, the mechanical parameters were determined as follows: the pressure strength was 12.578 MPa, the tensile strength was 1.274 MPa, the elastic modulus was 17.513GPa, and the Poisson's ratio was 0.219. To record the strain change near the fracture tip during the tests, strain gauges were pasted on both sides of the initial fracture after curing (Ni et al. 2015;Shi et al. 2020;Liu et al. 2020a, b;Zhang et al. 2020a, b, c, d).

Test equipment
According to Liu et al. (2018) and Shi et al. (2020), the three-field coupling triaxial test system was used to load water pressure and transverse pressure, and collect monitoring data, this system includes an universal testing machine of HS-3001A and a control computer. Among them, the universal testing machine of HS-3001A was used to load water pressure and transverse pressure (the load range of water pressure was 0 ~ 3.0 MPa, the load range of transverse pressure was 0 ~ 1000kN). The control computer can realize real-time acquisition of monitoring data. The test equipment is shown in Fig. 1.

Test scheme
The impacts of the angle between the transverse pressure stress and the initial fracture, the transverse pressure stress, the length and aperture of the initial fracture on the critical water pressure stress of HF are large, 4 levels of the above factors were selected for orthogonal experiment design, and a total of 16 groups of tests were designed. Additionally, in order to eliminate the influences of non-uniformity of materials and stress concentration on the test results, each test was carried out 3 times and the average values of the results were calculated (Huang et al. 2020;Yin et al. 2020). A total of 48 tests were conducted. The test scheme and results are shown in Table 1.

Model test results and analysis
When HF occurred on the cement mortar specimens, a large amount of water seeped from the fracture and made a dull sound of failure, the water pressure on the fracture plane decreased rapidly, but not completely penetrated, and there was still a certain residual strength of the specimen. At this time, the prefabricated fracture tip initiated micro-fractures, the direction of the fracture initiation and propagation was approximately perpendicular to the initial pressure stress direction. Figure 2 shows the duration curves of the transverse pressure stress and the water pressure stress of the I-1 test group. It can be seen that the water pressure stress in the first test dropped sharply from 1.981 MPa in the 1257 s to 0.394 MPa in the 1258 s. Thus HF occurred in the 1257 s, and the critical water pressure stress of HF was 1.981 MPa. Similarly, HF occurred in the second test and the third test at 1482 s and 1207 s, respectively, and the critical water pressure stress of HF was 2.056 MPa and 2.284 MPa.
Under the conditions of the category I (the angle between the transverse pressure stress and the initial fracture is 0°) and category II (the angle is 90°), the relationship between the transverse pressure stress and the critical water pressure stress of HF is shown in Formula 1 and Formula 2, respectively.
In the Formulas: ξ is the critical water pressure stress of HF, F 1 is the transverse pressure stress, δ is the mean square error. When the angle between the transverse pressure stress and the initial fracture is 0°, the transverse pressure stress promotes HF, at this time, the critical water pressure stress decreases with the increase of the transverse pressure stress. When the angle is 90°, the transverse pressure stress inhibits HF, at this time, the critical water pressure stress increases with the increase of the transverse pressure stress.
(1) = −0.1143F 1 + 2.2065 2 = 0.095 (2) = 0.085F 1 + 1.9905 2 = 0.104  Under the conditions of the category III, the relationship between the initial fracture length and the critical water pressure stress of HF is shown in Formula 3. Under the conditions of the category IV, the relationship between the initial fracture aperture and the critical water pressure stress of HF is shown in Formula 4.
In the Formulas: L is the initial fracture length, d is the initial fracture aperture. It can be seen that the larger the length and aperture of the initial fracture, the easier the HF occurs. When other conditions remain unchanged, the critical water pressure stress of initial fracture length equal 400 mm is on average 5.71% smaller than that of the initial fracture length equals 200 mm, while the critical water pressure stress of the initial fracture aperture equals 7 mm is on average 19.62% smaller than that of the initial fracture aperture equals 3 mm.

Numerical simulation method
To mutual authenticate the model test results, numerical simulations of the I-1 test group were carried out based on ABAQUS software. The researches on HF based on ABAQUS generally include two steps: Geostatic and Soil. Among them, Geostatic is the step of ground stress equilibrium analysis, aiming to construct an equilibriumed initial ground stress field and remove the deformation caused by initial ground stress. Soil is a seepagestress coupling analysis step to simulate fracture propagation under water pressure Sun et al. 2020a, b;Ni et al. 2015;Zeng et al. 2019). Due to the small size and light weight of the test model, it can satisfy the accuracy requirements without ground stress equilibrium analysis, so the numerical simulations only carried out Soil analysis. The three-dimensional numerical model of the I-1 test group consists of 50,500 grids and 55,539 nodes, as shown in Fig. 3.
In the simulation process, water was injected into the model from the injection orifice, which was located at the midpoint of the initial fracture, and the total injection time was 3000 s. The first 600 s were the injection speed-up stage, and the injection peak value was 0.06m 3 /min. The analysis step duration was set to 300, the number of incremental steps was 10,000, the initial incremental step was 0.01, the minimum incremental step was 1 × 10 -8 , the maximum incremental step was 10 and a single incremental step allowed a change in pore pressure of 100 MPa, other parameters are shown in Table 2.

Numerical simulation results and analysis
(1) Simulation results of the fracture aperture.
The fracture model, half-size model and full-size model are defined firstly. The fracture model is the most central 500 mm × 200 mm × 5 mm part of the model shown in Fig. 3, the  Fig. 3, and the full-size model is the model shown in Fig. 3. Figure 4 shows some representative simulation results of the fracture aperture of the fracture model, half-size model and full-size model at different incremental steps. At the beginning of the water injection (initial state), it was assumed that the initial fracture had no thickness (fracture aperture was 0). The fracture aperture was 4.270 mm at the 11th incremental step, while reached the maximum value of 5.126 mm when HF appeared at the 22nd incremental step, but the fracture plane had not reached fully penetrated at this time, which was consistent with the model test results.
(2) Simulation results of the mechanical parameters.  Figure 5 shows the pore pressure simulation results of the half-size model at different incremental steps. The pore pressure was 0 at the beginning of the water injection, while 27.02 MPa at the 11th incremental step. The pore pressure reached the maximum value of 0.1083GPa when HF appeared at the 22nd incremental step. Figure 6 shows the strain simulation results of the half-size model at different incremental steps. The strain was 0 at the beginning of the water injection, while 2.190 × 10 -2 at the 11th incremental step. The strain reached the maximum value of 5.626 × 10 -2 when HF appeared at the 22nd incremental step. Stiffness degradation means that the displacement of the peak point increases with the increase of the number of cycles when the same peak load is maintained under cyclic and repeated loads. Figure 7 shows the stiffness degradation index simulation results of the half-size model at different incremental steps. There was no stiffness degradation at the beginning of the water injection, the stiffness degradation index significantly increased at the 11th incremental step, and reached the maximum value when HF appeared at the 22nd incremental step.
It can be seen that the simulation results of the fracture aperture, pore pressure, strain and stiffness degradation index of the fracture model, half-size model and full-size model were basically consistent with the model test results, and were also basically consistent with the results of HF mechanism of rock revealed by other similar researches (Liu et al.   Li et al. 2021), which showed that the numerical simulation method was reasonable and correct.

HF mechanism analysis of rock
Griffith brittleness calculation criterion believes that due to the existence of micro-fractures in actual materials, when the nominal stress is still very low, the local stress concentration has reached a high value, exceeding the critical stress for fracture instability growth, so that the fracture expands rapidly and leads to brittle fracture (Huang et al. 2017). According to the Griffith brittleness calculation criterion, the premise of HF of brittle object is that the internal micro-fractures expand to macro fractures due to stress concentration, and the mechanical condition of brittle fracture is F Griffith < 0 (Liu et al. 2022a, b;Hou et al. 2022;Liu et al. 2022a, b). Many studies showed that HF mechanism did not fully comply with the Griffith brittleness calculation criterion because the rock (including cement mortar specimen) was not a standard brittle object, For example, Wang et al. (2006) studied the HF mechanism of cement mortar thick-walled cylinder specimens under different stress states, it was considered that the stress concentration induced by the non-uniformity of materials such as cement and sand was an important reason for HF. Liu et al. (2021) carried out the discrete element study on the interaction mechanism between HF and weak planes and considered that HF was the process of generating new fractures by the driving force of water in the weak planes. The water flow in the fracture has dual mechanical effects on the fracture wall: one is the hydrostatic pressure perpendicular to the fracture wall, and the other is the drag force along the flow direction (Liu et al. 2018). HF of rock is a quasi-brittle failure, including static stage, micro fracture propagation stage and macro fracture formation stage. The indicators that distinguish these three stages are the formation of micro fractures and macro fractures. Under the action of hydrostatic pressure, the initial fracture forms a deterioration zone at the fracture tip; with the increase of water pressure, the deterioration area is expanding under the dual mechanical effect of water flow, and micro fractures formed and expanded. When the stress intensity factors at the fracture tip exceed the fracture toughness, the micro-fractures expand into macro fractures, and HF occurs Li et al. 2022). The HF process of rock is shown in Fig. 8.   Fig. 7 Simulation results of the stiffness degradation index 1 3 5 Numerical simulations of rock slope failure under high groundwater pressure

Overview of the slope
Letuan-Qingshiguan section of G205 in China is located in the hinterland of central Shandong Mountains. The terrain of this area is dominated by mountains with complex geological structure, strong surface cutting and frequent human engineering activities. A rock slope along Letuan-Qingshiguan section is about 80°, and the height is about 22 m. The rock is mainly limestone, and the bedding structure is inclined. The slope is divided into the upper part, middle part and lower part by two large horizontal shear joints. The upper part is about 5 m thick, and the rock is highly weathered, the joints and fissures are developed, basically no vegetation grows, and there are obvious water erosion traces. There are three dangerous rocks on the upper part, which are distributed at 17 ~ 20 m high. The central part is about 8 m thick and moderately weathered, with a platform protruding outwards about 30 ~ 50 cm in width. The horizontal joints are well developed and the vegetation is dense. The lower part is about 9 m thick with a dip angle of about 20° and tends to the southeast with slight weathering and basically, no vegetation grows. The slope morphology and distribution of the dangerous rocks are shown in Fig. 9. The three-dimensional elevation data of the rock slope is plotted by Sufer software, as shown in Fig. 10, and the characteristic data of the dangerous rocks are shown in Table 3.  Chen et al. (2003) and Li et al. (2020) proposed three load combinations when studying the stability of dangerous rock, namely natural condition, rainstorm condition and seismic condition. The factors to be considered in each condition are as follows: 1. Natural condition: gravity and fracture water pressure in natural state. 2. Rainstorm condition: gravity and fracture water pressure in rainstorm state. 3. Seismic condition: gravity, fracture water pressure in natural state and seismic forces.
The stability coefficients under the above three conditions were calculated by taking the WYT3 as an example, and the selected parameters are shown in Table 4.
Through calculation, the stability coefficient of the WYT3 under the natural condition was 1.25, namely, it was in a stable state; the stability coefficients under the rainstorm condition and seismic condition were 0.96 and 1.07, namely, they were in an unstable state.

Numerical simulations of rock slope failure
According to the geological survey results, there is a high-pressure water-bearing body and a natural fracture plane inside the slope. The high-pressure water-bearing body is about  1.35 m in vertical direction and 0.95 m in horizontal direction, and the maximum water pressure stress is about 0.58 MPa. The natural fracture plane is formed by the action of the high-pressure water-bearing body, and its current length is about 3.307 m. The method proposed in Sect. 3 was used to simulate the whole processes of the expansion, penetration and failure of the natural fracture plane of WYT3 under the action of the high-pressure water-bearing body. In the numerical modeling, the natural condition was set, and the maximum principal stress criterion was selected. The failure mode was displacement mode, the failure displacement was 0.001 m, the viscosity regularization coefficient was 1.0 × 10 -4 , the elastic modulus was 15.0GPa, the Poisson's ratio was 0.25, the void ratio was 0.1, the permeability coefficient was 1.0 × 10 -7 and the soil expansion angle was 0.1. Other parameters are shown in Table 4. The model grids were all set to quadrilateral grids, including 586 grids and 642 nodes, as shown in Fig. 11. The water pressure stress F w applied by the high-pressure water-bearing body in the vertical direction was set to 0.6 MPa, the applied length was 1.35 m and the length of the natural fracture plane was 3.31 m. The numerical models of the high-pressure water-bearing body and the natural fracture plane are shown in Fig. 12.
According to the numerical analysis results, there are a large number of micro-fractures penetration in the process of natural fracture plane propagation, but these micro fractures are not the main cause of slope failure. Therefore, this paper only studies the development of the natural fracture plane (Zhuang et al. 2022). The natural fracture plane increased by 2.137 m under the action of high groundwater pressure, penetrated in 2467 s and occurred slope failure in 2468 s. Before the failure, the expansion of the natural fracture plane experienced three stages: slow development stage, rapid development stage and penetration stage. Among them, 0 ~ 1784s was the slow development stage, the natural fracture plane increased by 0.658 m, with an average growth rate of 0.000369 m/s; 1785 s ~ 2276 s was Fig. 11 Grid division diagram the rapid development stage, the natural fracture plane increased by 0.706 m, and the average growth rate was 0.001438 m/s; 2277s ~ 2467s was the penetration stage, the natural fracture plane increased by 0.773 m, and the average growth rate was 0.004068 m/s. The development of the natural fracture planes is shown in Fig. 13.

Conclusions
(1) The cement mortar specimens were selected to carry out HF model tests. To study the impacts of factors such as the angle between the transverse pressure stress and the initial fracture, the transverse pressure stress, the length and aperture of the initial fracture on the critical water pressure stress, 4 levels of the above factors were selected for orthogonal experiment design, and a total of 16 groups of tests were designed. Each test was carried out 3 times and a total of 48 tests were carried out. The results showed that when HF occurred on the cement mortar specimens, a large amount of water seeped from the fracture and made a dull sound of failure, the water pressure on the fracture plane decreased rapidly, but not completely penetrated, and there was still a certain residual strength of the specimen; The larger the length and aperture of the initial fracture, the easier the HF occurs; when the angle between the transverse pressure stress and the initial fracture is 0°, the transverse pressure stress promotes HF, while the angle is 90°, the transverse pressure stress inhibits HF.
(2) To mutual authenticate the model test results, numerical simulations of the I-1 test group were carried out, and analyzed the fracture aperture and mechanical parameters (including pore pressure, strain and stiffness degradation index) of the fracture model, half-size model and full-size model under different incremental steps, obtained that Natural fracture plane F w the simulation results were basically consistent with the model test results, and were also basically consistent with the results of HF mechanism of rock revealed by other similar researches, which showed that the numerical simulation method was reasonable and correct. (3) HF mechanism does not fully comply with the Griffith brittleness calculation criterion because the rock (including cement mortar specimen) is not a standard brittle object. The water flow in the fracture has dual mechanical effects on the fracture wall: one is the hydrostatic pressure perpendicular to the fracture wall, and the other is the drag force along the flow direction. HF of rock is a quasi-brittle failure, including static stage, micro fracture propagation stage and macro fracture formation stage. (4) There are 3 dangerous rocks in the rock slope along Letuotuan-Qingshiguan section of G205 in China, among which WYT3 is in a stable state under the natural condition and in an unstable state under the rainstorm condition and seismic condition. The fracture propagation of WYT3 under the action of high groundwater pressure experienced slow development stage, rapid development stage and penetration stage, among which the slow development stage lasted the longest, and the slope failure occurred immediately after the fracture penetrated (i.e., rock landslide occurred).