Research on the impact of pre‐existing geological fractures on hydraulic fracturing in high in situ stress environments

The formation of complex fracture network is a difficult problem in the process of deep rock reservoir reconstruction, and the key to the generation of intricate fracture systems in deep reservoirs is to communicate more pre‐existing fractures through hydraulic fracturing, so as to enhance oil recovery. At present, there are few studies on the interactions of hydraulic fractures with pre‐existing fractures under high stress. Therefore, this research studied the influence of the number of pre‐existing fractures on the characteristics of hydraulic fractures by using advanced hydraulic fracturing instruments, and then analyzed the characteristics of hydraulic fractures under different pre‐existing fracture number based on three‐dimensional topography scanning technology and three‐dimensional square box fractal dimension calculation method. Based on the three‐dimensional finite element discrete element method, further research is carried out. The laboratory test findings indicated that the increase of pre‐existing fractures will make the pump pressure curve more stable during the test, and will significantly improve the roughness of hydraulic fractures. The numerical simulation indicated that as the quantity of pre‐existing fractures grows, the fracture area and width also increase, and in instances of a substantial quantity of pre‐existing fractures, the hydraulic fractures are prone to bifurcation, which contributes to the development of intricate fracture networks. With the increase of natural fracture angle, the fracture width decreased, but the fracture area increased. In this study, the influence of the number of pre‐existing fractures on hydraulic fracturing was studied through laboratory tests and numerical simulations, which provided theoretical reference for engineering practice.

As science and technology progress, societal need for energy is on the rise, and shallow oil and gas resources are far from meeting the needs of production. 1 Therefore, exploring and developing deep unconventional reservoirs holds significant importance.The exploration and development of deep reservoirs is difficult because of the high in situ stress and the difficulty of fracture propagation. 2To form complex fracture networks in deep reservoirs and improve oil and gas recovery efficiency, making good use of pre-existing fractures in deep reservoirs is an effective way. 3 Hence, it is crucial to study the interaction between pre-existing fractures and hydraulic fractures.
Four methods are widely used in the study of reservoir rock hydraulic fracturing, which are theoretical research method, field test method, numerical simulation method, and laboratory test method. 4Theoretical research methods were popular when hydraulic fracturing was just coming out, but most of the theoretical research focused on shallow reservoirs, less on deep reservoirs, and many judgment criteria are not applicable to deep reservoirs. 5Field testing is the most accurate way to get an accurate understanding, but field testing is too expensive to conduct large-scale studies, 6 It may be used when studying issues related to geological hazards, such as the effects of hydraulic fracturing on faults. 7Hui et al. 8 studied the effects of hydraulic fracturing on seismic activity by collecting data from onsite hydraulic fracturing to develop fracture propagation models for unconventional hydraulic fracturing.Among all the methods, the most widely used method is the numerical simulation method, because it is simple, fast, low cost, and has low requirements for equipment and working environment, and considers many influencing factors, and can complete some operations that are difficult to achieve in conventional tests. 9The fourth method is the indoor test method, which is a similar test method.Through the indoor test method, the sample size and construction parameters are modified in proportion to the field test, which can obtain the test results similar to the field test, and the indoor test can obtain a complete fracture plane after hydraulic fracturing, which can facilitate a more in-depth exploration of the features of hydraulic fractures. 10he research of indoor hydraulic fracturing test mainly focus on the shallow reservoir rock, and deep reservoir rock was less studied. 11In addition, the samplemaking method is usually relatively simple, 12 and most of the hydraulic fracturing samples are cut and then spliced, which cannot properly restore the mechanical characteristics of hydraulic fractures in the formation, and the interface characteristics after splicing will greatly interfere with the expansion path of hydraulic fractures. 13In addition, even if the pre-existing fracture is embedded in the hydraulic fracturing sample, it is only a piece of paper or other material to replace the natural fracture, but this can only simulate the influence of microfractures, and can not study other complex preexisting fractures, let alone study the impact of preexisting fracture properties on hydraulic fracturing. 14here are four types of numerical simulation methods in hydraulic fracturing: extended finite element (XFEM), boundary element (BEM), discrete element (DEM), and finite element discrete element method (FDEM).XFEM characterizes discontinuity by adding a local enrichment function to the element containing the fracture, allowing the fracture to penetrate the mesh without the need for mesh reconstruction, so it can be widely used in hydraulic fracturing simulation tests, but it is not dominant in dealing with complex fracture networks. 15,16The BEM discretizes the domain's BEMs and approximates the boundary conditions using a function that fulfills the governing equation.However, the BEM shows a large discreteity in describing the propagation path of hydraulic fracture. 17,18DEM divides the research object into multiple units, 19 links each unit together through contact, and displays the overall mechanical characteristics by analyzing the contact between each unit. 20Therefore, DEM method has drawbacks in describing the action of fluid on solid. 21,22][25] The FDEM can not only reflect the advantages of high computational efficiency of the finite element method but also study the heterogeneous characteristics of rocks like the DEM.According to the fracture criteria, the joint element is fractured to simulate the initiation and propagation process of the fracture.Therefore, it is a good choice to use FDEM to study the fracture propagation characteristics during hydraulic fracturing.
In this paper, the influence of pre-existing fracture number on hydraulic fracturing is studied through laboratory experiments, and the fracture characteristics after hydraulic fracturing are studied by using the method for calculating the fractal dimension of a threedimensional (3D) cubic enclosure, and the influence mechanism of pre-existing fracture number on hydraulic fracture characteristics is analyzed.The influence of the number of Pre-Existing Fractures on the shape, area, and width of hydraulic fractures is further studied by the FDEM.In this paper, advanced test methods were used to study hydraulic fracturing, which can accurately reconstruct the mechanical environment of deep rock.The method of prefabrication of natural fractures was designed to study the influence of the number of natural fractures on hydraulic fracturing.A high-precision topography scanning system was used to obtain the digital information of hydraulic fracture surface, and advanced quantitative methods were used to evaluate the characteristics of hydraulic fracture.Theory, experiment, and simulation were combined to provide reference for engineering.

| TEST PROCEDURE 2.1 | Sample preparation
First, make natural fractures, wrap white ash with a paper bag, and then sprinkle water on the surface of the paper bag to make natural fractures.The hydraulic fracturing sample was poured with 42.5 R Portland cement, and different numbers of natural fractures were arranged during the pouring process, as depicted in Figure 1.A blind hole with a depth of 16.5 cm was drilled on the center surface of the specimen, and a slot with a width of 0.5 mm was cut at the bottom end of the blind hole with a depth of 15.5 cm as the initial fracture.

| Design of test program
Comparative tests SN17, SN20, SN21, and SN22 were conducted to study the effect of the number of fractures on hydraulic fracturing when the pre-existing fracture angle was 45°at σ V -σ H -σ h : 60-65-50 MPa.The injection rate Q is 20 mL/min. 26Fracturing fluid is slippery water, and the viscosity μ is 3 mPa•s. 27I G U R E 1 Sample preparation process.NF, natural fracture.

| Topography scanning of hydraulic fractures
The 3D scanner as depicted in Figure 2A was used to obtain the 3D fracture topography characteristics.

| Fractal dimension calculation of hydraulic fracture
Rock fissures in nature are complex structures.If Euclidean geometry is used to describe the points, lines, planes, volumes, and other features of rock joints, it is often difficult to use too many parameters.The fractal dimension derived from various box counting techniques can be represented using the following equation 29 : where D is the fractal dimension.δ i is the ith time square box count size.N is the number of square boxes (note: the fractal dimension is a double logarithmic fitting slope).
The measurement method of square box is used to calculate the fractal dimension, and the calculation principle is depicted in Figure 3.
In fact, when the fractal dimension of a complex geometry in a digital image is measured by the box dimension method, the smallest external box containing the complex geometry can be obtained first.In 3D space, the smallest external box can be taken as a cube or a cuboid.In a two-dimensional (2D) plane, the smallest external box is a square or a rectangle.By dividing the smallest external boxes equally, the number of boxes needed to cover the complex geometry with different sizes can be further obtained.Suppose the length, width, and height of the smallest external box are Δx, Δy, and Δz, respectively.If the smallest external box is equally divided, then the size δ i of the box covering the complex geometry can be expressed as In Formula (2), if Δx, Δy, and Δ z are not exactly equal, the smallest external box is a cuboid.If Δx = Δy = Δz, then the smallest external box is a cube.At the same time, if the complex geometry is on a 2D plane, then F I G U R E 4 Calculation method of hydraulic fracture area in three-dimensional space.

| Calculation method of hydraulic fracture area
Only the fractal dimension of hydraulic fracture surface can not describe the complexity of hydraulic fracture, so the hydraulic fracture area calculation method is introduced.First, the 3D scan data is converted into stl binary files.In the stl file, the 3D fracture surface is composed of several tiny triangular elements.In the stl file, the 3D coordinates of each vertex of each element triangle can be directly obtained, and then the side length of each triangular element can be obtained.Therefore, the area of each triangular element can be calculated one by one through Helen's formula, and the total area of the entire scanned fracture surface can be obtained by summing each triangular element, as shown in Figure 4.

| Analysis of test phenomenon
Evident from Figure 5, SN20, SN17, SN21, and SN22 samples have different numbers of natural fractures, but they all have natural fractures of the same angle (45°) made by 10 × 10 cm paper bag wrapped in white lime, and the σ V -σ H -σ h is 60-65-50 MPa, respectively.The pumping rate was 20 mL/min, and the fracturing fluid type was slippery water with viscosity μ = 3 mPa•s.Figure 5 is the profile comparison diagram.The natural fractures have been opened obviously, so they are A-Type crossing mode and secondary fractures appear.When the number of pre-existing fractures increases to 4, the way of hydraulic fractures penetrating pre-existing fractures is F-Type.When the number of hydraulic fractures increases to 6, the surface roughness of hydraulic fractures increases significantly.When there are eight pre-existing fractures, the hydraulic fractures penetrate the pre-existing fractures in a complex way.There are Type A and Type D near the wellbore.After Type D appears, there is an F-Type crossing mode at the edge of the sample.Therefore, when the triaxial stress is σ V -σ H -σ h : 60-65-50 MPa, with the increase of the number of natural fractures, the roughness of the fracture surface will also increase, and the mode of interactions between hydraulic fractures and natural fractures will become more complex, but the difference is that when the triaxial stress is σ V -σ H -σ h : 60-65-50 MPa, the roughness of the fracture surface is significantly higher, and the Type-D crossing mode will appear, and the fractures pass through both ends of the natural fractures, which promotes the development of intricate fracture networks.
Figure 6 shows the comparison of pumping pressure curves.When the triaxial stress is σ V -σ H -σ h : 60-65-50 MPa, the pumping pressure curve is significantly stable, especially in the process of forming a stable flow channel, and the pumping pressure curve is almost parallel to the X-axis.
Type-A in Figure 6 is the interaction pattern between hydraulic fracture and pre-existing fracture, natural fractures were cut through (natural fractures are not opened).And Type-F is that the hydraulic fracture passes through the pre-existing fracture at a certain point upward from end A. The Type-D mentioned earlier is a hydraulic fracture that passes through both ends of the natural fracture at the same time (the natural fracture opens).

| Fractal dimension analysis
Evident from Figure 7 (D stands for fractal dimension of hydraulic fracture), when the in situ stress is σ V -σ H -σ h : 60-65-50 MPa, the fractal dimension of the hydraulic fracture surface also grew in tandem with the rising count of pre-existing fractures.This indicated that the roughness of hydraulic fractures increased with the increase of the number of natural fractures.This is because, with the increase in the number of natural fractures, the tensile strength and shear strength of natural fractures and rock matrix composites are reduced, so that tensile fractures and shear fractures are easy to form under high stress and high pump pressure, and the roughness of the boundary between natural fractures and rock matrix is inherently high, so that the roughness of the fracture surface is improved.

| Fracture area analysis
Before calculating the hydraulic fracture area, it is first necessary to dig out the hydraulic fracture in the geomagic control, as shown in Figure 8.Then, the fracture area is calculated according to the method in 2.3.3.
It can be seen from Figure 8 that with the increase in the number of natural fractures, the area of natural fractures also increases, indicating that the increase of natural fractures can improve the reconstruction effect of the reservoir.

| DISCUSSIONS AND EXPANSIONS 4.1 | FDEM method based on cohesive
Due to the limitation of the study of fracture characteristics by laboratory test methods, 30 the influence of the number of Pre-Existing Fractures on the shape, area, and width of hydraulic fractures is further studied by FDEM.In the FDEM, the deformation behavior of rock element follows linear elastic criterion, and the seepage behavior of rock element follows Darcy law.The numerical solution of FDEM model is completed by using ABAQUS software. 31

| Fluid flow simulation equation
During fracturing, the flow flows from the opened cohesive units to the unopened cohesive units.The flow rate per unit length along the tangential direction is 32 where w is the cohesive unit thickness, m. µ is the fluid viscosity, mPa•s.p is the fluid pressure in the cohesive unit, MPa.| 941 The fluid transfer from the fracture to the adjacent reservoir rock can be calculated by 33 q c p p = ( − ), where q l is the local fluid loss per unit fracture surface area in the rock formation, m 3 .p m is the pore pressure in the formation, MPa.c l is the fluid loss coefficient, m/Pa•s.Darcy's law is adopted for fluid flow in rock mass 34 : where v is the average seepage velocity, m/s.Q is the flow rate, m 3 /s.A is the water crossing area, m 2 .k is the permeability coefficient, m/s.J is the hydraulic gradient.

| Damage model of interface unit
Before the damage occurs, the deformation behavior of cohesive interface unit is controlled by the law of traction and separation, as shown in Figure 9, satisfying the linear elastic relationship 35 : where ε δ T ε δ T ε δ T = / , = / , = / n s t n 0 s 0 t 0 .ε n is the normal strain.ε s , ε t is the strain of two transverse shear directions.T 0 is the initial thickness of the interface element, m.Δ n is the normal opening displacement, m. δ s , δ t is the opening displacement in both transverse shear directions, m.
When the maximum stress criterion is adopted, it is assumed that the cohesive unit breaks once the stress exceeds the critical threshold in any direction: where σ n 0 is the critical stress of cohesive element normal, MPa.σ σ , s t 0 0 is the critical stress of two tangents, MPa.<> indicates that cohesive element does not fail under pure compressive stress.
Taking bilinear constitutive as an example, the damage parameter changes of interface elements are shown in Figure 10.To ascertain the emergence of new fractures, one can monitor variations in the damage factor and the formula is as follows: Note that when the cohesive unit is compressed, where σ n is the predicted stress of cohesive element normal under the current strain according to the elastic tractor-separation criterion of undamaged line, MPa.σ σ ̅ , ̅ s t is the shear stress predicted by elastic tractorseparation criterion of undamaged line between two cohesive elements under the current strain, MPa.σ n , σ s , σ t corresponding to the actual stress in three directions, MPa.
When mixed-mode deformation occurs, the equivalent displacement can be expressed as When the nonlinear displacement expansion criterion is adopted, the damage factor in the stiffness degradation criterion is where β is the exponent describing the rate of fracture propagation.
F I G U R E 9 Tractive force and separation displacement relationship.

| Continuity equation and stress balance equation
Following the law of mass conservation, the influx of fluid into the rock mass over a specific time period equals the volume of fluid added to the rock mass.Therefore, the continuity equation is where K 0 is the result of multiplying the tensor representing the initial permeability coefficient with the water density.K r is the osmotic coefficient.K w is the bulk modulus, MPa.g is the acceleration of gravity, m/s 2 .n is the porosity, %.The stress-seepage coupling equation is Combined with the preceding formula, the calculation equation can be expressed as In this paper, the final simulation result is obtained by Newton iteration, and its time step will automatically increase or decrease according to convergence criteria.

| Model parameters
The corresponding 3D fracturing model is constructed to further study the effect of the number of pre-existing fractures on hydraulic fracture propagation during fracturing.The basic requirements of the model are as follows 31 : 1.The model size is 50 m × 50 m × 5 m. 2. The number of pre-existing fractures is the main variable in this simulation study.3. The model's outer boundary is characterized by fixed displacement and impermeability.4. The injection point and initial break unit are set in the center.5.The number of natural cracks was changed by changing the spacing between natural cracks.

| The effect of the number of natural fractures
In this section, discrete fracture networks are created by changing tangential and normal spacing between discrete fractures without increasing the number of bedding groups, so as to construct reservoir fracturing simulation models under the influence of different discrete fracture numbers and carry out numerical simulation analysis.In this simulation, the setting of tangential and normal spacing of pre-existing fractures is as follows: When the number of fractures is 82, the tangential spacing becomes between 2 and 3 m.When the number of fractures is 98, the tangential spacing becomes between 1 and 2 m.When the number of fractures is 137, the tangential spacing is between 1 and 2 m and the normal spacing is between 2 and 3 m.The corresponding results are shown in Figure 11.
Figure 11 shows the change of parameters under the influence of the number of pre-existing fractures.It can be seen from the above figure that the number of pre-existing fractures has a significant influence on fracture propagation under 3D conditions.With the increase in the number of pre-existing fractures, the time of pressure fractures reaching the boundary of the reservoir study area may advance or lag.At the same time, the corresponding phenomenon also appears in the fracture area, fracture volume, maximum fracture width, injection point fracture width, tensile failure ratio, and other parameters.It can be seen from Figure 11A that with the increase in the number of natural fractures, there is little influence on the initiation pressure.It can be seen from Figure 11B that with the increase in the number of pre-existing fractures, the fracture area gradually increases, which is consistent with the results of laboratory tests.It can be seen from Figure 11C,D that as the number of natural fractures increases, so does the width of the fractures.However, it can be found from Figure 11E,F that the relationship between the change in the number of natural fractures and the volume of fractures and the ratio of tensile failure is discrete.Therefore, further analysis should be carried out in combination with the results of fracture morphological differences, and the corresponding results are shown in Figure 12.It can be seen from Figure 12 that when the number of natural fractures is minimum, the phenomenon of artificial fractures spreading along the natural fractures occurs.However, this phenomenon is not significant.With the increase in the number of natural fractures, the intersection phenomenon of artificial fractures and natural fractures is more significant.Therefore, in the subsequent model with the largest number of pre-existing fractures, more obvious pre-existing fracture branching occurs.

| Influence of pre-existing fracture inclination
In this section, the corresponding influencing factors are analyzed by adjusting the inclination parameters in the discrete fracture network.The pre-existing fracture dip angles were set at 15°, 30°, 45°, 60°, and 75°, and the other parameters remained unchanged.The evolution of pressurized fracture parameters under the influence of pre-existing fracture inclination is shown in Figure 13.
Evident from Figure 13, for the 3D fracturing model, the pressure fracture appears to have reached the reservoir boundary more easily when the bedding inclination is 30°a nd 45°.With the increase of bedding dip angle, the difference in fluid pressure decreases.For fractured reservoirs with bedding angles of 15°, 60°, and 75°, With the increase of bedding dip angle, there is no significant law in the total fracture area when the reservoir pressure fracture reaches the reservoir boundary.The fracture width basically decreases with the increase of dip angle.By observing the variation of tensile failure ratio, it can be seen that the above phenomenon may be caused by the fact that reservoir pressure fractures are inclined to extend along pre-existing fractures when the inclination of pressure fractures increases, so the tensile failure ratio is largest when the inclination is 75°.Evident from Figure 13C,D with the increase of pre-existing fracture inclination, the fracture width decreases, because a lower fracture inclination will make the rock more prone to crack, increasing fracture width.However, it can be seen  from Figure 13E that a larger fracture inclination will increase the area of hydraulic fractures.This phenomenon can be explained by the criterion of energy release rate.Without considering the filtration, the energy dissipation effect of fluid is reflected in the increase of the width and length of hydraulic fractures, and the decrease of the fracture width naturally leads to the increase of the length, so the fracture area increases.
Figure 14 shows the difference in fracture morphology under the influence of fracture inclination.Evident from the above figure, pressure fractures in the 3D model tend to spread along pre-existing fractures.If and only if the pre-existing fracture inclination is large, the pressure fracture begins to show a pattern that tends to extend in the direction of the maximum principal stress.It can also be seen from Figure 14 that with the increase of fracture inclination, the fracture becomes more tortuous.It can also be seen from the figure that with the increase of fracture inclination, the fracture becomes more tortuous.

| CONCLUSIONS
In this study, the influence of the number of pre-existing fractures on hydraulic fracturing under high stress was studied by using laboratory test and numerical simulation methods, the characteristics of fractures were T A B L E 1 Fracturing model test scheme of samples containing pre-existing fractures. 28rial number Note: Table 2 shows the main parameters used in simulation models.
F I G U R E 14 The difference of fracture morphology under the influence of fracture inclination.
studied by 3D topography scanning technology combined with fractal theory, and a series of extended studies were conducted by using FDEM.The following conclusions were reached: 1. High-stress true triaxial hydraulic fracturing tests were conducted to study the effect of the number of preexisting fractures on hydraulic fracturing when the preexisting fracture angle was 45°.The number of fractures will affect the stability of pump pressure.The greater the number of fractures, the more stable the pump pressure and the easier it is to form a stable seepage channel.The fracture number will affect the bending degree of the fracture plane.The higher the number, the less bending degree the fracture plane will be affected, and the more stable the fracture expansion direction will be.With the increase in the number of pre-existing fractures, the roughness of fracture surface increased, and the interaction of hydraulic fractures with pre-existing fractures became more complicated.2. The numerical information of hydraulic fracture was obtained by 3D topography scanner, the fractal dimension of fracture surface was calculated by 3D square box fractal dimension calculation method, and the fracture area was calculated using the numerical information of the hydraulic fracture surface.The findings indicated that the increase of pre-existing fracture will significantly improve the roughness and the fracture area of hydraulic fracture, and the fractal dimension near the wellbore was significantly higher than that at other locations.3. The numerical simulation results of 3D FDEM showed that with the increase of the number of pre-existing fractures, the fracture area and width also increased, and when the number of pre-existing fractures was large, the intersection phenomenon of hydraulic fractures and preexisting fractures was more significant, and the hydraulic fractures were prone to bifurcation, which fostered the development of intricate fracture networks.With the increase of natural fracture angle, the fracture width decreased, but with the increase of natural fracture angle, the fracture area increased.4. The research showed that the existence of natural fractures in high-stress environment will help hydraulic fracturing to better transform the reservoir.Rock compressibility in areas with dense natural fractures is better, and the rational use of natural fractures is conducive to the formation of complex fracture networks.
The acquisition of fracture morphology.(A) Three-dimensional topography scanner.(B) Flow diagram.F I G U R E 3 Different square boxes covering three-dimensional fracture surface diagram.

F I G U R E 7
Curve of pre-existing fracture number-D.F I G U R E 8 Fracture area extraction map.CHANG ET AL.

F
I G U R E 10 Schematic diagram of damage parameter change of interface unit.

F
I G U R E 11 Parameter variation under the influence of the number of pre-existing fractures.(A) Time-fluid pressure, (B) time-fracture area, (C) time-injection point crack width, (D) time-maximum crack width, (E) fracture volume, and (F) time-tensile failure.

F
I G U R E 12The difference of fracture morphology under the influence of pre-existing fracture quantity.

F
I G U R E 13 Parameter variation under the influence of fracture inclination.(A) Time-fluid pressure, (B) time-fracture area, (C) time-injection point crack width, (D) time-maximum crack width, (E) fracture volume, and (F) time-tensile failure.
31in parameters used in simulation models.31 Note: Table1shows the specific scheme of the test design.Abbreviation: NF, natural fracture.TA B L E 2