Numerical investigation into the effect of natural fracture density on hydraulic fracture network propagation

Hydraulic fracturing is an important method to enhance permeability in oil and gas exploitation projects and weaken hard roofs of coal seams to reduce dynamic disasters, for example, rock burst. It is necessary to fully understand the mechanism of the initiation, propagation, and coalescence of hydraulic fracture network (HFN) caused by fluid flow in rock formations. In this study, a coupled hydro-mechanical model was built based on synthetic rock mass (SRM) method to investigate the effects of natural fracture (NF) density on HFN propagation. Firstly, the geometrical structures of NF obtained from borehole images at the field scale were applied to the model. Secondly, the micro-parameters of the proposed model were validated against the interaction between NF and hydraulic fracture (HF) in physical experiments. Finally, a series of numerical simulations were performed to study the mechanism of HFN propagation. In addition, confining pressure ratio (CPR) and injection rate were also taken into consideration. The results suggested that the increase of NF density drives the growth of stimulated reservoir volume (SRV), concentration area of injection pressure (CAIP), and the number of cracks caused by NF. The number of tensile cracks caused by rock matrix decrease gradually with the increase of NF density, and the number of shear cracks caused by rock matrix are almost immune to the change of NF density. The propagation orientation of HFN and the breakdown pressure in rock formations are mainly controlled by CPR. Different injection rates would result in a relatively big difference in the gradient of injection pressure, but this difference would be gradually narrowed with the increase of NF density. Natural fracture density is the key factor that influences the percentages of different crack types in HFN, regardless of the value of CPR and injection rate. The proposed model may help predict HFN propagation and optimize fracturing treatment designs in fractured rock formations. Disciplines Engineering | Science and Technology Studies Publication Details Chong, Z., Li, X., Chen, X., Zhang, J. & Lu, J. (2017). Numerical investigation into the effect of natural fracture density on hydraulic fracture network propagation. Energies, 10 (7), 914-1-914-33. This journal article is available at Research Online: http://ro.uow.edu.au/eispapers1/500


Introduction
Hydraulic fracturing is an important method used to enhance oil and gas production in horizontal wells [1,2].It can form hydraulic fracture networks (HFNs), along with natural fractures (NFs), in rock formations to areas far away from horizontal wells.The purpose of hydraulic fracturing is to stimulate the permeability of rock formations, and allow oil and gas to permeate into horizontal wells [3,4].In addition, a hard roof often leads to the danger of a hanging roof in the coal mining process, which, in turn, may result in dynamic disasters, such as rock burst [5,6].The most important solution is using hydraulic fracturing to weaken a hard roof.Therefore, it is important to investigate the effects of NF on HFN propagation.
Natural fractures are a kind of pre-existing flaw in rock formations, which determine the propagation trajectory of HFN in hydraulic fracturing (HF).The interaction between NF and HF has been widely investigated by many scholars [7][8][9].Renshaw and Pollard [10] proposed a simple criterion to identify crack propagation orientation during the mutually perpendicular interaction between NF and HF, and verified this criterion through physical experiments.Gu and Weng [11] further developed Renshaw's criterion and used the updated criterion to analyze non-orthogonal intersections between NF and HF.
The interaction between NFs and HFs in hydraulic fracturing can be divided into following six scenarios (Figure 1a,b) [10,12,13]: (1) Hydraulic fracturing directly crosses a NF (Figure 1c).Hydraulic fracturing can cross NFs easily when NFs have relatively strong cohesion, or the rock formations bear high confining stress.
In this case, the interaction interface opens with tension and the HF propagates along the same orientation with strong stress on the HF tip.(2) Hydraulic fracturing intersecting a NF (Figure 1d).When a HF crosses with a NF, fluid flows into the NF, but the HF does not initiate on the interface because of a non-existing weak face.
With the increase of fluid, the HF crosses the NF and propagates along the same orientation as before.As a result, they form T-shaped branches.(3) Hydraulic fracturing crossing with an offset (Figure 1e).When a HF interacts with a NF, a small offset occurs due to the stress concentration caused by localized interface separation and shear slip.(4) Hydraulic fracturing arrested by a NF (Figure 1f).When the strength of a NF is far smaller than that of the rock matrix, and the formations bear little confining pressure, the HF is arrested by the NF and cannot sufficiently transmit the NF because the interface fails in shear and slips.(5) Hydraulic fracturing shear slipping along a NF (Figure 1g).When the fluid pressure in a HF is smaller than the closure stress of a NF, the interface would be separated because of shear stress.The fluid permeability in the NF would be enlarged, leading to shear cracks propagating along the interface.(6) Hydraulic fracturing branching or turning at end of an NF (Figure 1h).When fluid pressure exceeds the closure stress of the interface, the NF would be widened by excessive tensile stress and allow fluid to flow in.Consequently, the NF would become a branch of the HFN.In addition, the concentration of fluid pressure would result in another HF initiation at the end of the NF.
Energies 2017, 10, 914 2 of 33 solution is using hydraulic fracturing to weaken a hard roof.Therefore, it is important to investigate the effects of NF on HFN propagation.Natural fractures are a kind of pre-existing flaw in rock formations, which determine the propagation trajectory of HFN in hydraulic fracturing (HF).The interaction between NF and HF has been widely investigated by many scholars [7][8][9].Renshaw and Pollard [10] proposed a simple criterion to identify crack propagation orientation during the mutually perpendicular interaction between NF and HF, and verified this criterion through physical experiments.Gu and Weng [11] further developed Renshaw's criterion and used the updated criterion to analyze non-orthogonal intersections between NF and HF.
The interaction between NFs and HFs in hydraulic fracturing can be divided into following six scenarios (Figure 1a,b) [10,12,13]: (1) Hydraulic fracturing directly crosses a NF (Figure 1c).Hydraulic fracturing can cross NFs easily when NFs have relatively strong cohesion, or the rock formations bear high confining stress.In this case, the interaction interface opens with tension and the HF propagates along the same orientation with strong stress on the HF tip.(2) Hydraulic fracturing intersecting a NF (Figure 1d).When a HF crosses with a NF, fluid flows into the NF, but the HF does not initiate on the interface because of a non-existing weak face.
With the increase of fluid, the HF crosses the NF and propagates along the same orientation as before.As a result, they form T-shaped branches.(3) Hydraulic fracturing crossing with an offset (Figure 1e).When a HF interacts with a NF, a small offset occurs due to the stress concentration caused by localized interface separation and shear slip.(4) Hydraulic fracturing arrested by a NF (Figure 1f).When the strength of a NF is far smaller than that of the rock matrix, and the formations bear little confining pressure, the HF is arrested by the NF and cannot sufficiently transmit the NF because the interface fails in shear and slips.(5) Hydraulic fracturing shear slipping along a NF (Figure 1g).When the fluid pressure in a HF is smaller than the closure stress of a NF, the interface would be separated because of shear stress.The fluid permeability in the NF would be enlarged, leading to shear cracks propagating along the interface.(6) Hydraulic fracturing branching or turning at end of an NF (Figure 1h).When fluid pressure exceeds the closure stress of the interface, the NF would be widened by excessive tensile stress and allow fluid to flow in.Consequently, the NF would become a branch of the HFN.In addition, the concentration of fluid pressure would result in another HF initiation at the end of the NF.Energies 2017, 10, 914 3 of 33 All the aforementioned scenarios assumed that the formations only have one NF.They are, however, abundant in real hydraulic fracturing process.In order to deal with this weakness, scholars proposed a discrete fracture network (DFN) model [14][15][16][17][18].A discrete fracture network is a network consisting of abundant discrete NFs, which distribute certain criteria.Liu et al. [19] evaluated nine factors affecting the equivalent permeability in a two-dimensional DFN model based on the mathematical formula.Huang et al. [20], who established a three-dimensional DFN model and compared it with the two-dimensional DFN model, also studied on factors affecting their equality.Meyer and Bazan [21] used an elliptical region containing two groups of parallel and uniformly spaced vertical fractures along directions of horizontal principal stress to represent the HFN in hydraulic fracturing.Weng et al. [13,22] built a special fracture model and studied the interaction criterion between NFs and HFs.
With these models, they successfully simulated the interactions between NFs and HFs.In addition, other numerical methods were also adopted to study HFNs in rock formations, including the boundary element method (BEM) [23,24], extended finite element method (XFEM) [25,26], and rock failure process analysis (RFPA) [27,28].The synthetic rock mass (SRM) method is one of many discrete element methods (DEMs), which is widely used to simulate the crack initiation, propagation, and coalescence in rock formations [29,30].In this model, rock mass is represented by an assembly of particles bound together, and pre-existing fractures are represented by DFN.In this model, the spatial properties of the fractures, obtained from geological statistics, are used as the input parameters of DFN.Although the models mentioned above can be used to study fractured rock formations, there are still many difficulties and challenges-the NF density influences the propagation of HFN in hydraulic fracturing is still unclear.Therefore, it is necessary to build a coupled hydro-mechanical model in fractured rock formations to investigate this issue.
In this paper, a SRM-based fracturing model is proposed to study the propagation of HFN in fractured rock formations.First, the geometrical structure of HFs used in this model was obtained from the statistical results of borehole images.Second, the micro-parameters of the model were calibrated through physical experiments that studied the interaction between NF and HF.The breakdown pressure of the model under different confining pressures was also compared with theoretical equations.Finally, a set of numerical simulations were performed to investigate the effects of NF density on HFN propagation.Additionally, external factors like confining pressure ratio (CPR) and injection rate were also considered.

Synthetic Rock Mass (SRM) Method
SRM method is considered as a proper model to study natural fracture (NF) in particle assemblies [31].Although SRM is a discontinuum method, particles can still bound together with the parallel bond method to simulate rock deformations.Tensile or shear cracks can initiate when the bonds between particles are broken.The movement and deformation of particles can influence the macro mechanical behavior of the model.Figure 2 shows contacts between particles, of which the force and force moment can be divided into two components [32]: where F i , F n , and F s are the total, normal-, and shear-directed force exerted by the parallel bond, respectively; M i , M n , and M s are the total, normal-, and shear-directed moment exerted by the parallel bond, respectively; and n i and t i are normal-and shear-directed unit vectors, respectively.

Figure 2.
The parallel bond and its breaks of inter-particle bond in SRM [32].
Friction coefficient pb μ of the contact between two particles is defined as: The increments of contact force and force moment are represented as: D are the increments of rotational angle which are added to the current value.
The tensile and shear strengths on particle contacts can be calculated with the following equations: where max σ and max τ are the tensile and shear strength, respectively; J is polar moment of inertia of the parallel bond cross section, respectively.SRM method is built by inserting discrete fracture network (DFN) into the adopted particle assembly.In the model, the NF is represented by smooth joint (SJ) logic (Figure 3).A SJ is used in SRM to simulate the crack development in rock mass [33].Once the joint plane is defined, the SJ would be distributed on the contacts between particles (Figure 3a).Each joint plane consists of two overlapping surfaces-Surface A and Surface B-particles can only move along these two surfaces Friction coefficient µ pb of the contact between two particles is defined as: The increments of contact force and force moment are represented as: where ∆F n and ∆F s are the increments of normal-and shear-directed force, respectively; k n and k s are the contact normal and shear stiffness, respectively; ∆U n and ∆U s are the increments of elastic force, which are added to current value; A and I are the area and moment of inertia of parallel bond cross section; ∆M n and ∆M s are the increments of normal-and shear-directed moment, respectively; and ∆θ n and ∆θ s are the increments of rotational angle which are added to the current value.
The tensile and shear strengths on particle contacts can be calculated with the following equations: where σ max and τ max are the tensile and shear strength, respectively; J is polar moment of inertia of the parallel bond cross section, respectively.SRM method is built by inserting discrete fracture network (DFN) into the adopted particle assembly.In the model, the NF is represented by smooth joint (SJ) logic (Figure 3).A SJ is used in SRM to simulate the crack development in rock mass [33].Once the joint plane is defined, the SJ would be distributed on the contacts between particles (Figure 3a).Each joint plane consists of two overlapping surfaces-Surface A and Surface B-particles can only move along these two surfaces (Figure 3b).The orientations of these surfaces are represented by joint unit normal vector nj which is defined through joint dip angle θ p .nj = (sin θ p , cos θ p ) (10) Energies 2017, 10, 914 5 of 33 (Figure 3b).The orientations of these surfaces are represented by joint unit normal vector ˆj n which is defined through joint dip angle p θ .The dot product of ˆj n and contact unit normal vector ˆc n are applied to determine the surface each particle lies.If

ˆ(sin ,cos )
, the particle lies in the nearer surface, such as Ball A lies in Surface A (Figure 3a).The existing bonds are removed to build SJ between particles.They are just like a set of elastic springs evenly distributed among the cross sections among particles.The area of SJ cross section A′ is defined as: The average radius R is defined as: min( , ) The force vector F and displacement vector U of SJ on the local coordinates of the joint plane are defined as: where n F and s F are normal force and shear force vectors of SJ, respectively, and n U and s U are normal relative displacement and shear relative displacement vectors of SJ.
The mechanical behavior of SJ follows the Coulomb sliding model, and the relative displacement increment of two particles interacting with each other in each time step is divided into two parts: normal and shear component.The normal and shear force of the SJ are updated using the following equations (Figure 4): where nj k and sj k are the normal and shear stiffness of SJ, respectively.Δ is an increment of the updated value.Maximum shear force s F * is determined by friction coefficient j μ of the SJ, defined as:  The dot product of nj and contact unit normal vector nc are applied to determine the surface each particle lies.If nj • nc > 0, the particle lies in the nearer surface, such as Ball A lies in Surface A (Figure 3a).The existing bonds are removed to build SJ between particles.They are just like a set of elastic springs evenly distributed among the cross sections among particles.The area of SJ cross section A is defined as: The average radius R is defined as: The force vector F and displacement vector U of SJ on the local coordinates of the joint plane are defined as: where F n and F s are normal force and shear force vectors of SJ, respectively, and U n and U s are normal relative displacement and shear relative displacement vectors of SJ.The mechanical behavior of SJ follows the Coulomb sliding model, and the relative displacement increment of two particles interacting with each other in each time step is divided into two parts: normal and shear component.The normal and shear force of the SJ are updated using the following equations (Figure 4): where k nj and k sj are the normal and shear stiffness of SJ, respectively.∆ is an increment of the updated value.Maximum shear force F * s is determined by friction coefficient µ j of the SJ, defined as: 4b).

Hydro-Mechanical Coupling in SRM
In SRM, the polygons formed through the contacts between particles are called fluid domains.Viscous fluid flow is based on the assumption that fluid domains are connected with fluid flow channels which refers to a series of contacts in particle assembly.Each fluid domain has its own volume which can be calculated by blue surrounding lines linking particle centers (Figure 5a).The fluid flow is driven by the pressure formed at the two ends of a fluid flow channel (Figure 5c-e).The fluid flow is controlled by Poiseuille equation which assumes that the fluid flow is laminar within the channel.Therefore, the volumetric flow rate Q can be defined as: where m is the fluid dynamic viscosity; P Δ is fluid pressure difference between the two neighboring domains; L and a are the length and aperture of the flow channel, respectively.

Hydro-Mechanical Coupling in SRM
In SRM, the polygons formed through the contacts between particles are called fluid domains.Viscous fluid flow is based on the assumption that fluid domains are connected with fluid flow channels which refers to a series of contacts in particle assembly.Each fluid domain has its own volume which can be calculated by blue surrounding lines linking particle centers (Figure 5a).The fluid flow is driven by the pressure formed at the two ends of a fluid flow channel (Figure 5c-e).The fluid flow is controlled by Poiseuille equation which assumes that the fluid flow is laminar within the channel.Therefore, the volumetric flow rate Q can be defined as: where µ is the fluid dynamic viscosity; ∆P is fluid pressure difference between the two neighboring domains; L and a are the length and aperture of the flow channel, respectively.

Hydro-Mechanical Coupling in SRM
In SRM, the polygons formed through the contacts between particles are called fluid domains.Viscous fluid flow is based on the assumption that fluid domains are connected with fluid flow channels which refers to a series of contacts in particle assembly.Each fluid domain has its own volume which can be calculated by blue surrounding lines linking particle centers (Figure 5a).The fluid flow is driven by the pressure formed at the two ends of a fluid flow channel (Figure 5c-e).The fluid flow is controlled by Poiseuille equation which assumes that the fluid flow is laminar within the channel.Therefore, the volumetric flow rate Q can be defined as: where m is the fluid dynamic viscosity; P Δ is fluid pressure difference between the two neighboring domains; L and a are the length and aperture of the flow channel, respectively.Figure 5b demonstrates that the fluid pressure (yellow spots) formed in the rock mass and the size of the spot is positively correlated to the magnitude of the pressure.In the calculation of fluid Figure 5b demonstrates that the fluid pressure (yellow spots) formed in the rock mass and the size of the spot is positively correlated to the magnitude of the pressure.In the calculation of fluid flow, the change of fluid domain volume was not taken into consideration in previous works to simplify hydro-mechanical coupling [34,35].The change, however, was considered in this paper.The increment of fluid pressure ∆P in the fluid domain can be calculated with the following equation: where K f is the fluid bulk modulus; ∆t is the increments of time step; V d is the pore volume; ∆V d is the change of value in the pore volume.Deformation of fluid domain is caused by the particle contact pressure, which in turn, leads to the deformation of whole rock mass.Consequently, the deformation of rock mass would drive the change of the contact force between particles.Cracks can initiate once contact stress exceeds the maximum tensile or shear stress in Equation (8) or Equation (9).In the proposed model, the aperture of the fluid channel can be calculated with an equation related to normal stress.An empirical equation has been adopted to describe the relationship between laminated pipe aperture and normal stress [36].
where a inf and a zeo are the aperture value at infinite and zero normal stress, respectively; coefficient ξ is the speed of aperture decay with the increasing of σ n and usually equals −0.15 [37].The units of a inf and a zeo are 'm'; the units of σ n is 'MPa'; coefficient ξ does not have any units.When σ n approaches infinity, a zeo would be infinitely close to a inf .The values of a zeo and a inf are defined as Equation ( 21) [34]: where k and V are the permeability and the total volume of rock mass.Once a crack initiates, the aperture of the fluid flow channel tends to be infinite, which may cause the problem of computational instability.Simulating the fluid flow between domains, which causes initiation of cracks, is vital to hydraulic fracturing.The cracks caused by the bond breakage between two particles can lead to instantaneous fluid flow.Fluid pressures P f in the two domains, after bond breakage, are considered as equivalent to the average of the pressures before bond breakage (Figure 6) [34,35,38].The updated fluid pressure calculation equation is defined as: where P f A and P f B are fluid pressures, related to Domain A and Domain B in Figure 6, respectively.This equation is not applicable when there is a large difference between two domain volumes.Therefore, once cracks initiate, the fluid pressure in the two related domains can be calculated through the following equation [39]: where ϕ is the model porosity; The modeling procedures of a hydro-mechanically coupled system is demonstrated in Figure 7.The deformation of rock formations would result in the change of fluid domain volume.At the same time, the aperture of the fluid flow channel would also change based on the change of particle contact stress.Because the fluid pressure is exerted on particle contact in the domain, the rock deformation and particle contact force can change at every time increment.In addition, a HFN consists of four types of cracks-namely, tensile crack of rock matrix (TCRM), shear crack of rock matrix (SCRM), tensile crack caused by natural fracture (TCNF), and shear crack caused by natural fracture (SCNF)-each of these were identified by an independent sub-program in this study.

Natural Fracture Structure
In order to quantitatively analyze NF structure, a borehole televiewer was performed on 9103 track roadways that were newly excavated in Huahong coal mine, China.The observation sites are shown in Figure 8a; the roof and ribs of track roadways are shown in Figure 8b,c, respectively.Five boreholes were observed with the ZKXG30 borehole images instrument (Figure 8d) on site, with which three boreholes (Φ28 × 10,000 mm) are on the roof and the other two (Φ28 × 5000 mm) on the ribs.From the observation, abundant NFs were found in borehole images (Figure 9).The modeling procedures of a hydro-mechanically coupled system is demonstrated in Figure 7.The deformation of rock formations would result in the change of fluid domain volume.At the same time, the aperture of the fluid flow channel would also change based on the change of particle contact stress.Because the fluid pressure is exerted on particle contact in the domain, the rock deformation and particle contact force can change at every time increment.In addition, a HFN consists of four types of cracks-namely, tensile crack of rock matrix (TCRM), shear crack of rock matrix (SCRM), tensile crack caused by natural fracture (TCNF), and shear crack caused by natural fracture (SCNF)-each of these were identified by an independent sub-program in this study.The modeling procedures of a hydro-mechanically coupled system is demonstrated in Figure 7.The deformation of rock formations would result in the change of fluid domain volume.At the same time, the aperture of the fluid flow channel would also change based on the change of particle contact stress.Because the fluid pressure is exerted on particle contact in the domain, the rock deformation and particle contact force can change at every time increment.In addition, a HFN consists of four types of cracks-namely, tensile crack of rock matrix (TCRM), shear crack of rock matrix (SCRM), tensile crack caused by natural fracture (TCNF), and shear crack caused by natural fracture (SCNF)-each of these were identified by an independent sub-program in this study.

Natural Fracture Structure
In order to quantitatively analyze NF structure, a borehole televiewer was performed on 9103 track roadways that were newly excavated in Huahong coal mine, China.The observation sites are shown in Figure 8a; the roof and ribs of track roadways are shown in Figure 8b,c, respectively.Five boreholes were observed with the ZKXG30 borehole images instrument (Figure 8d) on site, with which three boreholes (Φ28 × 10,000 mm) are on the roof and the other two (Φ28 × 5000 mm) on the ribs.From the observation, abundant NFs were found in borehole images (Figure 9).

Natural Fracture Structure
In order to quantitatively analyze NF structure, a borehole televiewer was performed on 9103 track roadways that were newly excavated in Huahong coal mine, China.The observation sites are shown in Figure 8a; the roof and ribs of track roadways are shown in Figure 8b,c, respectively.Five boreholes were observed with the ZKXG30 borehole images instrument (Figure 8d) on site, with which three boreholes (Φ28 × 10,000 mm) are on the roof and the other two (Φ28 × 5000 mm) on the ribs.From the observation, abundant NFs were found in borehole images (Figure 9).Natural fracture plane (NFP) in a small area can be regarded as a two-dimension plane.An intersection line is formed while drilling the borehole on the plane.By unfolding the borehole wall, the intersection line presents three geometrical shapes.The vertical intersection between borehole and NFP forms one horizontal intersection line (Figure 10a), parallel intersection forms two vertical intersection lines (Figure 10b) and inclined intersection forms one standard sine curve (Figure 10c).
The dip angle of NFP α can be calculated by the following equation:  Natural fracture plane (NFP) in a small area can be regarded as a two-dimension plane.An intersection line is formed while drilling the borehole on the plane.By unfolding the borehole wall, the intersection line presents three geometrical shapes.The vertical intersection between borehole and NFP forms one horizontal intersection line (Figure 10a), parallel intersection forms two vertical intersection lines (Figure 10b) and inclined intersection forms one standard sine curve (Figure 10c).
The dip angle of NFP α can be calculated by the following equation: Natural fracture plane (NFP) in a small area can be regarded as a two-dimension plane.An intersection line is formed while drilling the borehole on the plane.By unfolding the borehole wall, the intersection line presents three geometrical shapes.The vertical intersection between borehole and NFP forms one horizontal intersection line (Figure 10a), parallel intersection forms two vertical Energies 2017, 10, 914 10 of 33 intersection lines (Figure 10b) and inclined intersection forms one standard sine curve (Figure 10c).The dip angle of NFP α can be calculated by the following equation: where h ab is the height difference of Points A and B in Figure 10c; r is the radius of the borehole; x A and x B are the x-coordinate of Points A and B, respectively.
Energies 2017, 10, 914 10 of 33 where ab h is the height difference of Points A and B in Figure 10c; r is the radius of the borehole;  Although the theoretical quantitative analysis is conducted on the structure of NFP obtained from borehole images in Figure 10, the uncertainties in the borehole imaging process make the analysis much more complicated in real situations.The impact of gravel stripping and collapse within borehole imaging is inevitable.However, considering NF density-the main factor considered in this paper is determined on the basis of the statistics of abundant NFs-minor errors have little impact on the statistical results of NFs [40].
In this paper, NF density was determined by borehole images of accumulative length of 40 m.One-thousand-and-seventy-four NFs, in total, were observed, and the geometrical parameters of each fracture were calculated, including length, aperture, and dip angle.Figure 11 shows the borehole images, three-dimensional virtual maps, and statistical results of natural fractures which were obtained from a borehole on the right side of the roof plate.
According to statistical results (Figure 12), the frequency is over 240 (the highest) from 180 mm to 200 mm, and all lengths of NFs were normally distributed.The equation is defined as: The R-square of Equation ( 25) is 0.93.Although NF aperture also varied to a certain degree, the variance was small compared to the variance of the length.Therefore, the average of the aperture values of 1074 NFs was selected as the parameter of the aperture parameter (d = 2.56 mm) for the model.The dip angles of NF were found to be normally distributed in three areas.Natural fracture density ρ can be defined as the sum of NF lengths sum l to the measured area sum A ratio: Although the theoretical quantitative analysis is conducted on the structure of NFP obtained from borehole images in Figure 10, the uncertainties in the borehole imaging process make the analysis much more complicated in real situations.The impact of gravel stripping and collapse within borehole imaging is inevitable.However, considering NF density-the main factor considered in this paper is determined on the basis of the statistics of abundant NFs-minor errors have little impact on the statistical results of NFs [40].
In this paper, NF density was determined by borehole images of accumulative length of 40 m.One-thousand-and-seventy-four NFs, in total, were observed, and the geometrical parameters of each fracture were calculated, including length, aperture, and dip angle.Figure 11 shows the borehole images, three-dimensional virtual maps, and statistical results of natural fractures which were obtained from a borehole on the right side of the roof plate.
According to statistical results (Figure 12), the frequency is over 240 (the highest) from 180 mm to 200 mm, and all lengths of NFs were normally distributed.The equation is defined as: The R-square of Equation ( 25) is 0.93.Although NF aperture also varied to a certain degree, the variance was small compared to the variance of the length.Therefore, the average of the aperture values of 1074 NFs was selected as the parameter of the aperture parameter (d = 2.56 mm) for the model.The dip angles of NF were found to be normally distributed in three areas.Natural fracture density ρ can be defined as the sum of NF lengths l sum to the measured area A sum ratio: Based on the statistical results, NF density in this paper was determined as 48/m (ρ = 48/m).The results were used to generate the NF in the proposed coupled hydro-mechanical model through an independent sub-program.
The results were used to generate the NF in the proposed coupled hydro-mechanical model through an independent sub-program.The results were used to generate the NF in the proposed coupled hydro-mechanical model through an independent sub-program.

Model Procedure
In this paper, a 1.0 m × 1.0 m model was built, in which natural fracture (NF) was generated according to the statistical results in Section 3. In order to research the effects of NF on HFN propagation, ±25% and ±50% changes were made to the statistical results of NF density during the process of generating NF using the model, which means that NF density ρ was set as 24 (−50%), 36 (−25%), 48, 60 (+25%) and 96/m (+50%), respectively.The length and dip angles of NFs, however, followed a normal distribution based on the statistical results of geometrical structure in Section 3. The average NF aperture value (d = 2.56 mm) remained the same.
The injection hole was formed through removing particles in the model center.The confining pressure in vertical orientation is represented by H S and horizontal orientation by h S .The stress difference of confining pressure h Δ ranges from 0 to 18 MPa.

Model Validation
In the proposed model, the micro-properties of intact specimens were firstly calibrated through a trial and error method.According to the authors' previous works [41], the calibrated microproperties are listed in Table 1.The comparison between the numerical simulation results and physical experimentation results of rock formation under different CPRs is shown in Figure 13.The figure shows that the results obtained from numerical simulation are in good agreement with those from the physical experiments.

Model Procedure
In this paper, a 1.0 m × 1.0 m model was built, in which natural fracture (NF) was generated according to the statistical results in Section 3. In order to research the effects of NF on HFN propagation, ±25% and ±50% changes were made to the statistical results of NF density during the process of generating NF using the model, which means that NF density ρ was set as 24 (−50%), 36 (−25%), 48, 60 (+25%) and 96/m (+50%), respectively.The length and dip angles of NFs, however, followed a normal distribution based on the statistical results of geometrical structure in Section 3. The average NF aperture value (d = 2.56 mm) remained the same.
The injection hole was formed through removing particles in the model center.The confining pressure in vertical orientation is represented by S H and horizontal orientation by S h .The stress difference of confining pressure ∆h ranges from 0 to 18 MPa.

Model Validation
In the proposed model, the micro-properties of intact specimens were firstly calibrated through a trial and error method.According to the authors' previous works [41], the calibrated micro-properties are listed in Table 1.The comparison between the numerical simulation results and physical experimentation results of rock formation under different CPRs is shown in Figure 13.The figure shows that the results obtained from numerical simulation are in good agreement with those from the physical experiments.
Table 1.Micro-properties used in intact rock for simulated specimens after calibration [41].

Micro-Parameters Unit Values
The    The micro-parameters of the SJ were calibrated based on the interactive behavior between NF and HF in physical experiments.The selected intersection angles θ and ∆h (∆h = |S H − S h |), in numerical simulations, were the same to those conducted in physical experiments [42][43][44].When the SJ and fluid parameters were fixed, θ and ∆h in the rock formation played a decisive role in the results of interaction between HF and NF.In this study, three typical θ (θ = 30 • ; 60 • ; 90 • ) values were adopted.The SJ and fluid parameters were adjusted repeatedly with the use of the trial-and-error method, so that the results of the interactions between HFs and NFs in the numerical simulations were the same as those in the physical experiments.Figure 14a,b,d,e,g,h demonstrates that the HF opens the NF; Figure 14c,f,i demonstrates that the HF crosses the NF.The physical experiments and numerical simulation suggested that: (1) The possibility of HFs crossing NFs increases as θ increases.This is due to the increasing difficulty for fluid to flow into NFs when they are oriented closer to the normal value of the HF propagation path.(2) The possibility of HFs crossing NFs increases along with the increase of ∆h.When there is a greater pressure difference between the horizontal orientation and vertical orientation, it is easier for HF propagation along the maximum principal stress.
The separate line curve presented a negative exponent change.The calibrated scenarios appeared similar to the experimental findings of [42][43][44], as shown in Figure 15.The SJ parameters after calibration and fluid parameters of numerical simulation are listed in Tables 2 and 3.
Energies 2017, 10, 914 14 of 33 The physical experiments and numerical simulation suggested that: (1) The possibility of HFs crossing NFs increases as θ increases.This is due to the increasing difficulty for fluid to flow into NFs when they are oriented closer to the normal value of the HF propagation path.
(2) The possibility of HFs crossing NFs increases along with the increase of h Δ .When there is a greater pressure difference between the horizontal orientation and vertical orientation, it is easier for HF propagation along the maximum principal stress.
The separate line curve presented a negative exponent change.The calibrated scenarios appeared similar to the experimental findings of [42][43][44], as shown in Figure 15.The SJ parameters after calibration and fluid parameters of numerical simulation are listed in Tables 2 and 3.In addition to the initiation, propagation, and coalescence of HFN, breakdown pressure was also taken into consideration in simulations of hydraulic fracture with the proposed model.A typical regression equation was proposed to calculate the breakdown pressure.Wang et al. [45] analyzed the effect of tensile strength and initial stress parameters on breakdown pressure.They discussed the coefficients of the equation, based on fitting curves, and suggested that both anisotropy and inhomogeneity of rocks would result in different equation coefficients.In addition, Fjar et al. [46] put Comparison between the interactive behavior in the physical experiment and numerical simulation.In addition to the initiation, propagation, and coalescence of HFN, breakdown pressure was also taken into consideration in simulations of hydraulic fracture with the proposed model.A typical regression equation was proposed to calculate the breakdown pressure.Wang et al. [45] analyzed the effect of tensile strength and initial stress parameters on breakdown pressure.They discussed the coefficients of the equation, based on fitting curves, and suggested that both anisotropy and inhomogeneity of rocks would result in different equation coefficients.In addition, Fjar et al. [46] put forward another regression equation in which the Poisson's ratio υ of the crack was taken into consideration to calculate breakdown pressure P b : where ω is the poroelastic coefficient; σ t is the rock tensile strength.
In order to verify the applicability of the proposed model, the simulated breakdown pressure value was compared with the calculated value through Equation ( 27) under same confining ratio in this study.The tensile strength of the NF density model built in Section 3 was 2.48 MPa (ρ = 48/m), υ was 0.22, and ω ranged from 0.15 to 0.24 [47].The initial pore pressure was set as 0 (P 0 = 0), the confining pressure on the y axis (S H ) was set as 15 MPa, and on the x-axis (S h ) varied from 10 to 15 MPa.The parameters of rock and fluid were adopted from Tables 1-3. Figure 16 demonstrates the comparison between the result obtained from theoretical equation and that from numerical simulations under different CPRs.The errors between the two were within 15%, meaning that the proposed model is applicable to hydraulic fracture simulation.
where ω is the poroelastic coefficient; t σ is the rock tensile strength.
In order to verify the applicability of the proposed model, the simulated breakdown pressure value was compared with the calculated value through Equation ( 27) under same confining ratio in this study.The tensile strength of the NF density model built in Section 3 was 2.48MPa ( ρ = 48/m), υ was 0.22, and ω ranged from 0.15 to 0.24 [47].The initial pore pressure was set as 0 ( 0 P = 0), the confining pressure on the y axis (SH) was set as 15 MPa, and on the x-axis (Sh) varied from 10 to 15 MPa.The parameters of rock and fluid were adopted from Tables 1-3. Figure 16 demonstrates the comparison between the result obtained from theoretical equation and that from numerical simulations under different CPRs.The errors between the two were within 15%, meaning that the proposed model is applicable to hydraulic fracture simulation.

Results and Analysis
In addition to NF density, some other external factors also influence the initiation, propagation, and coalescence of HFN:CPR and injection rate.Based on the calibrated micro-parameters in Section 4, a series of numerical simulations were performed to research the influence of NF density on the evolution of HFN propagation.Meanwhile, the influence of CPR and injection rate were also taken into consideration.

Evolution Process
Stimulated reservoir volume (SRV) is widely used to describe the efficiency of rock reservoir fracturing, which is defined as [48]: where x L and y L are the width and the length of SRV, respectively; the thickness of the numerical model z H is assumed to be that of the unit length.
Figure 17 demonstrates the evolution of SRV under different NF densities, showing that a constant declining of SRV increases rate.Under different NF density conditions, the increase of time step (>2000) would lead to bigger variances in SRV.The higher the NF density, the bigger the SRV.

Results and Analysis
In addition to NF density, some other external factors also influence the initiation, propagation, and coalescence of HFN:CPR and injection rate.Based on the calibrated micro-parameters in Section 4, a series of numerical simulations were performed to research the influence of NF density on the evolution of HFN propagation.Meanwhile, the influence of CPR and injection rate were also taken into consideration.

Evolution Process
Stimulated reservoir volume (SRV) is widely used to describe the efficiency of rock reservoir fracturing, which is defined as [48]: where L x and L y are the width and the length of SRV, respectively; the thickness of the numerical model H z is assumed to be that of the unit length.Figure 17 demonstrates the evolution of SRV under different NF densities, showing that a constant declining of SRV increases rate.Under different NF density conditions, the increase of time step (>2000) would lead to bigger variances in SRV.The higher the NF density, the bigger the SRV.Particularly, SRV under maximum NF density (ρ = 72/m) was 31.54%larger than that of the minimum NF density (ρ = 24/m).
Particularly, SRV under maximum NF density ( ρ = 72/m) was 31.54%larger than that of the minimum NF density ( ρ = 24/m).Figure 18 demonstrates the evolution of injection pressure under different NF densities.Within 800 time steps, the concentration area of injection pressure (CAIP) remained similar under different NF densities, because rock itself had resistance capacity to deformation in this stage.With the increase of time steps, however, the gradient of injection pressure and CAIP gradually differentiated from their prior values, because the strength of the rock with high NF density was low, and even low injection pressure could fracture the rock formations.The higher the NF density, the smaller of the injection pressure gradient, and the bigger the CAIP.In addition, high NF density made the propagation and coalescence of cracks easier in rock formations, which resulted in a larger CAIP. Figure 18 demonstrates the evolution of injection pressure under different NF densities.Within 800 time steps, the concentration area of injection pressure (CAIP) remained similar under different NF densities, because rock itself had resistance capacity to deformation in this stage.With the increase of time steps, however, the gradient of injection pressure and CAIP gradually differentiated from their prior values, because the strength of the rock with high NF density was low, and even low injection pressure could fracture the rock formations.The higher the NF density, the smaller of the injection pressure gradient, and the bigger the CAIP.In addition, high NF density made the propagation and coalescence of cracks easier in rock formations, which resulted in a larger CAIP.  Figure 18 demonstrates the evolution of injection pressure under different NF densities.Within 800 time steps, the concentration area of injection pressure (CAIP) remained similar under different NF densities, because rock itself had resistance capacity to deformation in this stage.With the increase of time steps, however, the gradient of injection pressure and CAIP gradually differentiated from their prior values, because the strength of the rock with high NF density was low, and even low injection pressure could fracture the rock formations.The higher the NF density, the smaller of the injection pressure gradient, and the bigger the CAIP.In addition, high NF density made the propagation and coalescence of cracks easier in rock formations, which resulted in a larger CAIP.

Considering Confining Pressure
In this section, the influence of NF density on HFN propagation under different CPRs, S H : S h = 10:15, 10:10, and 10:5 were investigated respectively, and all parameters required were adopted from Tables 1-3. Figure 20 demonstrates the influence of NF density on injection pressure under different CPRs.When S H : S h = 10:15, the injection pressure in the rock with the NF density of 24/m approached the maximum value (60 MPa), because both rock strength and confining pressure reached their maximum levels.All CAIPs of rock formations with different densities were mainly developed horizontally.As to S H : S h = 10:10 or 10:5, the injection pressure decreased sharply.When horizontal confining pressure and vertical confining pressure were same (S H : S h = 10:10), however, CAIP expanded gradually and distributed evenly on horizontal and vertical orientations.When confining pressure on vertical orientation was larger than that on horizontal orientation (S H : S h = 10:5), CAIP mainly developed along vertical orientation.Under the same CPR, although CAIP remained almost constant, the gradient of maximum injection pressure decreased gradually with the increase of NF density.Particularly, injection pressure presented the greatest decline during the process of NF density-increasing from ρ = 24/m to 48/m.The fitting strength envelope curves reflected the influence of NF density on breakdown pressure under different CPRs, and their corresponding equations are demonstrated in Figure 21 and Table 5, respectively.When NF density and vertical confining pressure remained constant, the increase of horizontal confining pressure would result in the growth of breakdown pressure.The maximum breakdown pressure approached 40 MPa when the horizontal confining pressure reached 15 MPa.The increase of NF density would narrow the variances in breakdown pressures.Under the same CPR, the breakdown pressure presented negative exponential distribution with an increase of NF density, which meant that there was little variance when NF density was high.

Considering Confining Pressure
In this section, the influence of NF density on HFN propagation under different CPRs, SH: Sh = 10:15, 10:10, and 10:5 were investigated respectively, and all parameters required were adopted from Tables 1-3. Figure 20 demonstrates the influence of NF density on injection pressure under different CPRs.When SH: Sh = 10:15, the injection pressure in the rock with the NF density of 24/m approached the maximum value (60 MPa), because both rock strength and confining pressure reached their maximum levels.All CAIPs of rock formations with different densities were mainly developed horizontally.As to SH: Sh = 10:10 or 10:5, the injection pressure decreased sharply.When horizontal confining pressure and vertical confining pressure were same (SH: Sh = 10:10), however, CAIP expanded gradually and distributed evenly on horizontal and vertical orientations.When confining pressure on vertical orientation was larger than that on horizontal orientation (SH: Sh = 10:5), CAIP mainly developed along vertical orientation.Under the same CPR, although CAIP remained almost constant, the gradient of maximum injection pressure decreased gradually with the increase of NF density.Particularly, injection pressure presented the greatest decline during the process of NF density-increasing from ρ = 24/m to 48/m.The fitting strength envelope curves reflected the influence of NF density on breakdown pressure under different CPRs, and their corresponding equations are demonstrated in Figure 21 and Table 5, respectively.When NF density and vertical confining pressure remained constant, the increase of horizontal confining pressure would result in the growth of breakdown pressure.The maximum breakdown pressure approached 40 MPa when the horizontal confining pressure reached 15 MPa.The increase of NF density would narrow the variances in breakdown pressures.Under the same CPR, the breakdown pressure presented negative exponential distribution with an increase of NF density, which meant that there was little variance when NF density was high.When NF density was relatively low (ρ = 24/m), it was TCRM that mainly distributed both at the border of SRV and around the injection hole.SCNF only initiated in concentration area of TCRM, and intersection or coalescence was found between the two.Under this condition, NF density was low and the strength of the rock formation itself was enough, which stopped HFN from propagating to areas far away from the injection hole.With the increase in NF density (ρ = 36/m or 48/m), TCRMs still dominated in number around the injection hole.Since the increase of NF density led to significant growth in the number of cracks caused by natural fracture in HFN, SCNF mainly distributed at the border of the SRV.When NF density reached a high level (ρ = 60/m or 72/m), SCNF became dominant, both at the border of the SRV and around the injection hole.TCRM was only sparsely distributed in HFN, but the number of these cracks was still far more than that of SCRM and TCNF.
Figure 23 demonstrates the incremental accumulation of TCRM, SCRM, TCNF, and SCNF under different CPRs during fracturing process.In the rock formations with the same NF densities, the total number of cracks under S H : S h = 10:10 was less than that in the other two CPRs.The results suggested that as the CPR approached 1.0, it became more difficult for the fracturing treatment to form rocks.When ρ = 24/m, TCRM had the most value under all three CPRs.With an increase of NF density (ρ = 36/m), SCNF values grew gradually and reached roughly an equal level with TCRM value.When ρ ≥ 48/m, the SCNF value overtook TCRM value and became dominant among different crack types.The higher the NF density, the bigger of the gap between the value of SCNF and those of the other crack types.The value of TCNF also grew gradually with the increase in NF density.When ρ = 72/m, the SCNF value almost doubled TCRM value.In particular, NF density and confining pressure had little influence on the SCRM value.When NF density was relatively low ( ρ = 24/m), it was TCRM that mainly distributed both at the border of SRV and around the injection hole.SCNF only initiated in concentration area of TCRM, and no intersection or coalescence was found between the two.Under this condition, NF density was low and the strength of the rock formation itself was enough, which stopped HFN from propagating to areas far away from the injection hole.With the increase in NF density ( ρ = 36/m or 48/m), TCRMs still dominated in number around the injection hole.Since the increase of NF density led to significant growth in the number of cracks caused by natural fracture in HFN, SCNF mainly distributed at the border of the SRV.When NF density reached a high level ( ρ = 60/m or 72/m), SCNF became dominant, both at the border of the SRV and around the injection hole.TCRM was only sparsely distributed in HFN, but the number of these cracks was still far more than that of SCRM and TCNF.
Figure 23 demonstrates the incremental accumulation of TCRM, SCRM, TCNF, and SCNF under different CPRs during fracturing process.In the rock formations with the same NF densities, the total number of cracks under SH: Sh = 10:10 was less than that in the other two CPRs.The results suggested that as the CPR approached 1.0, it became more difficult for the fracturing treatment to form rocks.The envelope curves reflect the relationship between NF density and the total number of cracks under different CPRs are shown in Figure 24 and Table 6. Figure 24 reveals an exponential relationship between the total number of cracks and NF density.When NF density was relatively low ( ρ = 24/m or 36/m), different CPRs resulted in large differences in the total number of cracks, which reached the maximum level when SH: Sh = 10:5.The differences in the total number of cracks under different CPRs were gradually narrowed with an increase of NF density ( ρ ≥ 36/m).In other words, the variance of CPR had little influence on the change in the total number of cracks in rock formations with a high NF density.

Considering Injection Rate
In this section, the influence of NF density on HFN propagation under different injection rates: 2.0 × 10 −6 m 3 /s, 2.5×10 −6 m 3 /s, and 3.0×10 −6 m 3 /s, was investigated, and all required parameters were The fitting envelope curves reflect the relationship between NF density and the total number of cracks under different CPRs are shown in Figure 24 and Table 6. Figure 24 reveals an exponential relationship between the total number of cracks and NF density.When NF density was relatively low (ρ = 24/m or 36/m), different CPRs resulted in large differences in the total number of cracks, which reached the maximum level when S H : S h = 10:5.The differences in the total number of cracks under different CPRs were gradually narrowed with an increase of NF density (ρ ≥ 36/m).In other words, the variance of CPR had little influence on the change in the total number of cracks in rock formations with a high NF density.The fitting envelope curves reflect the relationship between NF density and the total number of cracks under different CPRs are shown in Figure 24 and Table 6.

Considering Injection Rate
In this section, the influence of NF density on HFN propagation under different injection rates: 2.0 × 10 −6 m 3 /s, 2.5×10 −6 m 3 /s, and 3.0×10 −6 m 3 /s, was investigated, and all required parameters were

Considering Injection Rate
In this section, the influence of NF density on HFN propagation under different injection rates: 2.0 × 10 −6 m 3 /s, 2.5×10 −6 m 3 /s, and 3.0×10 −6 m 3 /s, was investigated, and all required parameters were adopted from Tables 1-3. Figure 25 demonstrates the influence of NF density on injection pressure under different injection rates.was minimum under the injection rate of 2.0 × 10 −6 m 3 /s, and with the increase of injection rate, the CAIP expanded gradually, especially along horizontal orientation.Meanwhile, the maximum injection pressure also grew with an injection rate increase.When the injection rate reached 3.0 × 10 −6 m 3 /s, the injection pressure in the rock with ρ = 24/m became the biggest among all specimens and the CAIP also expanded with an increase of NF density.Meanwhile, under an injection rate of 3.0 × 10 −6 m 3 /s, the expansions of CAIP in the horizontal and vertical orientation were almost the same as in the rock formation, with ρ = 72/m.The fitting strength envelope curves reflect the influence of NF density on breakdown pressure under different injection rates, and their corresponding equations are demonstrated in Figure 26 and Table 7, respectively.The results revealed a negative exponential relationship between breakdown pressure and NF density.When ρ = 24/m, different injection rates resulted in the biggest difference in breakdown pressure.Particularly, under the injection rate of 3.0 × 10 −6 m 3 /s, the breakdown pressure doubled under the injection rate of 2.0 × 10 −6 m 3 /s.The difference in breakdown pressure gradually became smaller with an increase of NF density.This difference almost disappeared when ρ = 72/m, which meant that the rock formation with this NF density was fragile.
Energies 2017, 10, 914 23 of 33 adopted from Tables 1-3. Figure 25 demonstrates the influence of NF density on injection pressure under different injection rates.CAIP was minimum under the injection rate of 2.0 × 10 −6 m 3 /s, and with the increase of injection rate, the CAIP expanded gradually, especially along horizontal orientation.Meanwhile, the maximum injection pressure also grew with an injection rate increase.
When the injection rate reached 3.0 × 10 −6 m 3 /s, the injection pressure in the rock with ρ = 24/m became the biggest among all specimens and the CAIP also expanded with an increase of NF density.
Meanwhile, under an injection rate of 3.0 × 10 −6 m 3 /s, the expansions of CAIP in the horizontal and vertical orientation were almost the same as in the rock formation, with ρ = 72/m.
The fitting strength envelope curves reflect the influence of NF density on breakdown pressure under different injection rates, and their corresponding equations are demonstrated in Figure 26 and Table 7, respectively.The results revealed a negative exponential relationship between breakdown pressure and NF density.When ρ = 24/m, different injection rates resulted in the biggest difference in breakdown pressure.Particularly, under the injection rate of 3.0 × 10 −6 m 3 /s, the breakdown pressure doubled under the injection rate of 2.0 × 10 −6 m 3 /s.The difference in breakdown pressure gradually became smaller with an increase of NF density.This difference almost disappeared when ρ = 72/m, which meant that the rock formation with this NF density was fragile.During the primary stage of fluid injection, all injection rates were high and exceeded the threshold of crack initiation.Therefore, although injection rates were different, the type, number, and distribution of cracks around injection holes on the rock formations were similar.With further propagation, however, the distribution of cracks under different injection rates tended to be different.The area of crack propagation under the minimum injection rate of 2.0 × 10 −6 m 3 /s also became the   During the primary stage of fluid injection, all injection rates were high and exceeded the threshold of crack initiation.Therefore, although injection rates were different, the type, number, and distribution of cracks around injection holes on the rock formations were similar.With further propagation, however, the distribution of cracks under different injection rates tended to be different.The area of crack propagation under the minimum injection rate of 2.0 × 10 −6 m 3 /s also became the Figure 27 reveals the influence of NF density on the crack distribution under different injection rates.During the primary stage of fluid injection, all injection rates were high and exceeded the threshold of crack initiation.Therefore, although injection rates were different, the type, number, and distribution of cracks around injection holes on the rock formations were similar.With further propagation, however, the distribution of cracks under different injection rates tended to be different.The area of crack propagation under the minimum injection rate of 2.0 × 10 −6 m 3 /s also became the smallest one.Nevertheless, when NF density was high enough (ρ ≥ 60/m), this phenomenon disappeared and crack propagation areas under all injection rates tended to be similar.Additionally, the injection rate had no influence on the percentage of different types of cracks.From high NF density to low NF density, the dominant values for crack types also shifted from TCRM to SCNF.
Energies 2017, 10, 914 25 of 33 smallest one.Nevertheless, when NF density was high enough ( ρ ≥ 60/m), this phenomenon disappeared and crack propagation areas under all injection rates tended to be similar.Additionally, the injection rate had no influence on the percentage of different types of cracks.From high NF density to low NF density, the dominant values for crack types also shifted from TCRM to SCNF.The fitting envelope curves reflect the relationship between NF density and total number of cracks under different injection rates are shown in Figure 29 and Table 8, respectively.The fitting envelope curves reflect the relationship between NF density and total number of cracks under different injection rates are shown in Figure 29 and Table 8, respectively.Figure 29 reveals an increase in the total number of cracks with higher injection rates and NF density.The increase in injection rates can result in a linear growth of the crack number in the fracturing process.reveals an increase in the total number of cracks with higher injection rates and NF density.The increase in injection rates can result in a linear growth of the crack number in the fracturing process.

Particle Size
In Section 5, the size of the model was 1.0 × 1.0 m, and the average particle size was 9 mm.During the NF process in fractured rock formation, the sensitivity of the particle size was analyzed to determine the influence on HFN.In the 1.0 × 1.0 m model, five particle sizes of 4.50 (−50%), 6.25 (−25%), 9.0, 11.25 (+25%) and 13.50 (+50%) mm were selected for NF examinations.
Figure 30 demonstrates the influence of different particle sizes on injection pressure and crack distribution of the model when the NF density is 48/m and SH: Sh = 10:5.When the model size remained unchanged, the variance in particle size had little influence on the injection pressure, while CAIP still developed along the maximum principal stress orientation.Additionally, the distribution of maximum injection pressures was essentially similar.Since the variance of particle size changed the number of contacts in the model, there was some difference in the quantity of cracks.Particle size had little influence, however, on the trend of fracture propagation and the proportion of crack types.
Ding [49] suggested that when the ratio of model size (L) to the average particle size (R) is small, the difference of particle size would lead to a bigger variance in macro-mechanical behavior of the model.When L/R ≥ 50, however, the variance of particle size has limited influence on the macro-mechanical behavior of the model.In this study, the average particle size was from 4.50 to 13.50 mm and L/R ratio was bigger than 50.These figures show that the particle size had a limited effect on mechanical behavior when the L/R ratio was greater than 50, which agrees well with published results [49].In Section 5, the size of the model was 1.0 × 1.0 m, and the average particle size was 9 mm.During the NF process in fractured rock formation, the sensitivity of the particle size was analyzed to determine the influence on HFN.In the 1.0 × 1.0 m model, five particle sizes of 4.50 (−50%), 6.25 (−25%), 9.0, 11.25 (+25%) and 13.50 (+50%) mm were selected for NF examinations.
Figure 30 demonstrates the influence of different particle sizes on injection pressure and crack distribution of the model when the NF density is 48/m and S H : S h = 10:5.When the model size remained unchanged, the variance in particle size had little influence on the injection pressure, while CAIP still developed along the maximum principal stress orientation.Additionally, the distribution of maximum injection pressures was essentially similar.Since the variance of particle size changed the number of contacts in the model, there was some difference in the quantity of cracks.Particle size had little influence, however, on the trend of fracture propagation and the proportion of crack types.
Ding [49] suggested that when the ratio of model size (L) to the average particle size (R) is small, the difference of particle size would lead to a bigger variance in macro-mechanical behavior of the model.When L/R ≥ 50, however, the variance of particle size has limited influence on the macro-mechanical behavior of the model.In this study, the average particle size was from 4.50 to 13.50 mm and L/R ratio was bigger than 50.These figures show that the particle size had a limited effect on mechanical behavior when the L/R ratio was greater than 50, which agrees well with published results [49].

Model Size
In order to study the influence of model size on the effect of fractures, six sizes of 1.0, 10, 25, 50, 75 and 100 m were selected.Figure 31 shows the influence of different model sizes on the injection pressure and crack distribution of the model when the NF density was 48/m and S H : S h = 10:5.The variance in model size had a remarkable effect on injection pressure and crack distribution.The injection pressure greatly increased with the increase in model size.When the model size reached 10 m, cracks appeared in the NF of rock formations.Thereafter, cracks gradually appeared in the rock matrix.Under the premise of different model sizes, however, CAIP and fractures still developed along the vertical orientation, i.e., the maximum principal stress orientation.
Based on many numerical simulation models, our explanation for the discrepancy among the results was that our model was larger in size than the models in Section 5, and that the size of the model can influence the value of the injection pressure and crack distribution, which is our next subject of research.
Another reason was the differentiation of parameters input in the model.In Section 4, the input parameters were calibrated in line with physical experiments at the laboratory scale.If these parameters were applied to a model of a larger size, size effects and anisotropy would inevitably appear, leading to different injection pressures and fracture distributions.Thus, it is necessary to study the size effects and anisotropy before the proposed model goes through HF under different sizes in order to determine REV.This will be our next topic of research.

Model Size
In order to study the influence of model size on the effect of fractures, six sizes of 1.0, 10, 25, 50, 75 and 100 m were selected.Figure 31 shows the influence of different model sizes on the injection pressure and crack distribution of the model when the NF density was 48/m and SH: Sh = 10:5.The variance in model size had a remarkable effect on injection pressure and crack distribution.The injection pressure greatly increased with the increase in model size.When the model size reached 10 m, cracks appeared in the NF of rock formations.Thereafter, cracks gradually appeared in the rock matrix.Under the premise of different model sizes, however, CAIP and fractures still developed along the vertical orientation, i.e., the maximum principal stress orientation.
Based on many numerical simulation models, our explanation for the discrepancy among the results was that our model was larger in size than the models in Section 5, and that the size of the model can influence the value of the injection pressure and crack distribution, which is our next subject of research.
Another reason was the differentiation of parameters input in the model.In Section 4, the input parameters were calibrated in line with physical experiments at the laboratory scale.If these parameters were applied to a model of a larger size, size effects and anisotropy would inevitably appear, leading to different injection pressures and fracture distributions.Thus, it is necessary to study the size effects and anisotropy before the proposed model goes through HF under different sizes in order to determine REV.This will be our next topic of research.

Conclusions
A series of numerical modellings for hydraulic fracturing were performed using a modified coupled hydro-mechanical model in SRM to investigate the effect of NF density on HFN propagation.In addition, two key factors-CPR and injection rate-were also analyzed and the following conclusions could be drawn: (1) Natural fracturing density had a significant influence on SRV, CAIP, and the numbers of different types of cracks in the HFN forming process.For a specimen with a high NF density, hydraulic fracturing often resulted in a larger SRV, larger CAIP, and more SCNF and TCNF.The number of TCRM, however, decreased gradually with an increase in NF density, and SCRM is almost immune to NF density.(2) In the fracturing process, it was the CPR that controlled the differences in HFN propagation orientation and breakdown pressure.With an increase in NF density, the differences gradually narrowed.When confining pressures were the same in the horizontal and vertical orientations, the number of cracks initiated in the HFN were minimal.(3) A difference in injection rate can lead to a greater variance in the gradient of injection pressure in rock formations.In addition, the change in injection rate can only influence the HFN in areas far away from the injection hole, and this influence gradually disappeared with an increase in NF density ( ρ ≥ 60/m).

Conclusions
A series of numerical modellings for hydraulic fracturing were performed using a modified coupled hydro-mechanical model in SRM to investigate the effect of NF density on HFN propagation.In addition, two key factors-CPR and injection rate-were also analyzed and the following conclusions could be drawn: (1) Natural fracturing density had a significant influence on SRV, CAIP, and the numbers of different types of cracks in the HFN forming process.For a specimen with a high NF density, hydraulic fracturing often resulted in a larger SRV, larger CAIP, and more SCNF and TCNF.The number of TCRM, however, decreased gradually with an increase in NF density, and SCRM is almost immune to NF density.(2) In the fracturing process, it was the CPR that controlled the differences in HFN propagation orientation and breakdown pressure.With an increase in NF density, the differences gradually narrowed.When confining pressures were the same in the horizontal and vertical orientations, the number of cracks initiated in the HFN were minimal.(3) A difference in injection rate can lead to a greater variance in the gradient of injection pressure in rock formations.In addition, the change in injection rate can only influence the HFN in areas far away from the injection hole, and this influence gradually disappeared with an increase in NF density (ρ ≥ 60/m).
(4) Regardless of the CPR and injection rate, the percentages of different crack types in HFN were principally controlled by NF density.In rock formations with a high enough NF density (ρ ≥ 72/m), the total number of cracks under different CPRs and injection rates were similar, which was conducive to parameter optimization in fracture treatment.
It is crucial for a reservoir stimulation to understand the effects of NF density on HFN propagation, as well as HF and NF interactions.Furthermore, the proposed model will cause scale effects and anisotropy, with an increase in size-this topic will be investigated in future studies.

Figure 1 .
Figure 1.Geometry of a hydraulic fracture network (HFN) depicting various scenarios of interaction with natural fractures (NFs).

Figure 1 .
Figure 1.Geometry of a hydraulic fracture network (HFN) depicting various scenarios of interaction with natural fractures (NFs).

Figure 2 .
Figure 2. The parallel bond and its breaks of inter-particle bond in SRM [32].

Figure 3 .
Figure 3. Smooth joint (SJ) contact model in two dimensions.

Figure 3 .
Figure 3. Smooth joint (SJ) contact model in two dimensions.

Figure 4 .
Figure 4. Smooth joint force-displacement law (a) normal force versus normal displacement; (b) shear force versus shear displacement.

Figure 5 .
Figure 5. Synthetic rock mass (SRM) method used in coupled hydro-mechanical simulation (a) fluid domains (blue polygons), centers of domains (red circles), and flow pipes (red lines) composing the fluid network; (b) distribution of fluid pressure; (c) domain and flow channel; (d) fluid pressure build up; (e) change of aperture.

Figure
Figure 5b demonstrates that the fluid pressure (yellow spots) formed in the rock mass and the size of the spot is positively correlated to the magnitude of the pressure.In the calculation of fluid

Figure 4 .
Figure 4. Smooth joint force-displacement law (a) normal force versus normal displacement; (b) shear force versus shear displacement.

Figure 4 .
Figure 4. Smooth joint force-displacement law (a) normal force versus normal displacement; (b) shear force versus shear displacement.

Figure 5 .
Figure 5. Synthetic rock mass (SRM) method used in coupled hydro-mechanical simulation (a) fluid domains (blue polygons), centers of domains (red circles), and flow pipes (red lines) composing the fluid network; (b) distribution of fluid pressure; (c) domain and flow channel; (d) fluid pressure build up; (e) change of aperture.

Figure 5 .
Figure 5. Synthetic rock mass (SRM) method used in coupled hydro-mechanical simulation (a) fluid domains (blue polygons), centers of domains (red circles), and flow pipes (red lines) composing the fluid network; (b) distribution of fluid pressure; (c) domain and flow channel; (d) fluid pressure build up; (e) change of aperture.

Figure 6 .
Figure 6.Bond breakage and fluid pressure balancing according to Equation (23) in the two domains: (a) fluid pressure before a crack initiation and (b) fluid pressure after a crack initiation.

Figure 7 .
Figure 7. Modeling procedures for a hydro-mechanical coupling.

Force
Modifying the fluid pipes; Modifying the volume and domain Fluid induced force calculation Failure criteria satisfied?Tensile crack of rock matrix (TCRM) Shear crack of rock matrix (SCRM) Tensile crack caused by natural fracture (TCNF) Shear crack caused by natural fracture (TCNF)

Figure 6 .
Figure 6.Bond breakage and fluid pressure balancing according to Equation (23) in the two domains: (a) fluid pressure before a crack initiation and (b) fluid pressure after a crack initiation.

Figure 6 .
Figure 6.Bond breakage and fluid pressure balancing according to Equation (23) in the two domains: (a) fluid pressure before a crack initiation and (b) fluid pressure after a crack initiation.

Figure 7 .
Figure 7. Modeling procedures for a hydro-mechanical coupling.

Figure 8 .
Figure 8.(a) Plan view of the study site at Huahong coal mine, China; photographs of (b) roof (c) rib; (d) test equipment of borehole image.

Figure 9 .
Figure 9. Optical televiewer image of the borehole drilled vertically in the roof at the Huahong coal mine, China.

Figure 8 . 33 Figure 8 .
Figure 8.(a) Plan view of the study site at Huahong coal mine, China; photographs of (b) roof (c) rib; (d) test equipment of borehole image.

Figure 9 .
Figure 9. Optical televiewer image of the borehole drilled vertically in the roof at the Huahong coal mine, China.

Figure 9 .
Figure 9. Optical televiewer image of the borehole drilled vertically in the roof at the Huahong coal mine, China.
Ax and Bx are the x-coordinate of Points A and B, respectively.

Figure 10 .
Figure 10.The intersection of the drill hole and the structure surface: (a) vertical intersection; (b) parallel intersection and (c) inclined intersection.

Figure 10 .
Figure 10.The intersection of the drill hole and the structure surface: (a) vertical intersection; (b) parallel intersection and (c) inclined intersection.

Figure 11 .
Figure 11.(a) Typical point; (b) virtual map of borehole image in 3-D; (c) plan view of the drill hole; (d) measurements of natural fracture (NF).

Figure 11 .
Figure 11.(a) Typical point; (b) virtual map of borehole image in 3-D; (c) plan view of the drill hole; (d) measurements of natural fracture (NF).

Figure 11 .
Figure 11.(a) Typical point; (b) virtual map of borehole image in 3-D; (c) plan view of the drill hole; (d) measurements of natural fracture (NF).
minimum particle radius (R min ) mm 1.0 Ratio of maximum and minimum particle radius (R max /R min ) modulus of the parallel bond (E c ) GPa 14.4 Normal stiffness of the parallel bond (mean) MPa 14.7 Normal stiffness of the parallel bond (std deviation) MPa 3.6 Shear stiffness of the parallel bond (mean) MPa 9.2 Shear stiffness of the parallel bond (std deviation) MPa 2.4

Figure 13 .
Figure 13.Comparison between experimental and numerical deviation stress versus axial strain curves.The micro-parameters of the SJ were calibrated based on the interactive behavior between NF and HF in physical experiments.The selected intersection angles θ and h Δ (

Figure 14 .
Figure 14.Numerical results of HF and NF interactions at different values of θ and h Δ .

Figure 13 .
Figure 13.Comparison between experimental and numerical deviation stress versus axial strain curves.

Energies 2017, 10 , 914 13 of 33 Figure 13 .
Figure 13.Comparison between experimental and numerical deviation stress versus axial strain curves.The micro-parameters of the SJ were calibrated based on the interactive behavior between NF and HF in physical experiments.The selected intersection angles θ and h Δ (

Figure 14 .
Figure 14.Numerical results of HF and NF interactions at different values of θ and h Δ .

Figure 14 .
Figure 14.Numerical results of HF and NF interactions at different values of θ and ∆h.

Figure 15 .Table 2 .Table 3 .
Figure 15.Comparison between the interactive behavior in the physical experiment and numerical simulation.Table2.Micro-properties used in smooth joint (SJ) contact for simulated specimens after calibration.ParametersUnit Values

Figure 15 .
Figure 15.Comparison between the interactive behavior in the physical experiment and numerical simulation.

Figure 16 .
Figure 16.Comparison of breakdown pressure under different confining pressure ratios obtained by the theoretical equation and numerical simulations.

Figure 16 .
Figure 16.Comparison of breakdown pressure under different confining pressure ratios obtained by the theoretical equation and numerical simulations.

Figure 17 .
Figure 17.Change of stimulated reservoir volume (SRV) with time steps when natural fracture (NF) density ranges from 24 to 72/m when injection rate is 2.0 × 10 −6 m 3 /s and S H : S h = 10:5.

Figure 18 .
Figure 18.Evolution process of injection pressure induced by NF densities of rock formations when injection rate is 2.0 × 10 −6 m 3 /s and SH: Sh = 10:5.

Figure 19
Figure19shows the changes of different crack types under five NF densities.The fitting curve equations reflecting the relationship between the number of cracks and NF density are listed in Table4.When NF density was relatively low ( ρ = 24/m), TCRM had the highest value under all confining ratios.The medium NF density ( ρ = 36/m or 48/m) would lead to the increase of TCNF and SCNF values, and SCNF value presented an exponential increase in trend.In this process, TCRM and SCNF were close to each other in value.The value of SCNF, however, would exceed that of TCRM and become dominant when NF density was high ( ρ = 60/m or 72/m).Additionally, NF density had limited influence on the value of SCRM, which always remained minimal compared to other crack types.

Figure 19 .
Figure 19.Variation of tensile crack of rock matrix (TCRM), shear crack of rock matrix (SCRM), tensile crack of natural fracture (TCNF), and shear crack of rock matrix (SCNF) in various NF densities.

Table 4 .Figure 18 .
Figure 18.Evolution process of injection pressure induced by NF densities of rock formations when injection rate is 2.0 × 10 −6 m 3 /s and S H : S h = 10:5.

Figure 19 33 Figure 18 .
Figure19shows the changes of different crack types under five NF densities.The fitting curve equations reflecting the relationship between the number of cracks and NF density are listed in Table4.When NF density was relatively low (ρ = 24/m), TCRM had the highest value under all confining ratios.The medium NF density (ρ = 36/m or 48/m) would lead to the increase of TCNF and SCNF values, and SCNF value presented an exponential increase in trend.In this process, TCRM and SCNF were close to each other in value.The value of SCNF, however, would exceed that of TCRM and become dominant when NF density was high (ρ = 60/m or 72/m).Additionally, NF density had limited influence on the value of SCRM, which always remained minimal compared to other crack types.

Figure 19
Figure19shows the changes of different crack types under five NF densities.The fitting curve equations reflecting the relationship between the number of cracks and NF density are listed in Table4.When NF density was relatively low ( ρ = 24/m), TCRM had the highest value under all confining ratios.The medium NF density ( ρ = 36/m or 48/m) would lead to the increase of TCNF and SCNF values, and SCNF value presented an exponential increase in trend.In this process, TCRM and SCNF were close to each other in value.The value of SCNF, however, would exceed that of TCRM and become dominant when NF density was high ( ρ = 60/m or 72/m).Additionally, NF density had limited influence on the value of SCRM, which always remained minimal compared to other crack types.

Figure 19 .
Figure 19.Variation of tensile crack of rock matrix (TCRM), shear crack of rock matrix (SCRM), tensile crack of natural fracture (TCNF), and shear crack of rock matrix (SCNF) in various NF densities.

Table 4 .
Fitting curve equations of the number of cracks for different NF densities.Confining Pressure Ratios Fitting Curve Equations R-Squared Tensile crack of rock matrix (TCRM) 0.18 32.2 y x = − + 0.96 Shear crack of rock matrix (SCRM) 0.67 y = 0.98 Tensile crack caused by natural fracture (TCNF) 0.06 1.5 y x = + 0.81 Shear crack caused by natural fracture (SCNF)

Figure 19 .
Figure 19.Variation of tensile crack of rock matrix (TCRM), shear crack of rock matrix (SCRM), tensile crack of natural fracture (TCNF), and shear crack of rock matrix (SCNF) in various NF densities.

Figure 20 .
Figure 20.Change of injection pressure with CPR when NF density ranges from 24/m to 72/m and injection rate is 2.0 × 10 −6 m 3 /s.

Figure 21 .Table 5 .
Figure 21.Variation of breakdown pressure with NF density under different confining pressure ratios and injection rate is 2.0 × 10 −6 m 3 /s.

Figure 20 . 33 Figure 20 .
Figure 20.Change of injection pressure with CPR when NF density ranges from 24/m to 72/m and injection rate is 2.0 × 10 −6 m 3 /s.

Figure 21 .Table 5 .
Figure 21.Variation of breakdown pressure with NF density under different confining pressure ratios and injection rate is 2.0 × 10 −6 m 3 /s.

Figure 22
Figure22demonstrates the influence of NF density on crack distribution under different CPRs.When NF density was relatively low (ρ = 24/m), it was TCRM that mainly distributed both at the border of SRV and around the injection hole.SCNF only initiated in concentration area of TCRM, and intersection or coalescence was found between the two.Under this condition, NF density was low and the strength of the rock formation itself was enough, which stopped HFN from propagating to areas far away from the injection hole.With the increase in NF density (ρ = 36/m or 48/m), TCRMs still dominated in number around the injection hole.Since the increase of NF density led to significant growth in the number of cracks caused by natural fracture in HFN, SCNF mainly distributed at the border of the SRV.When NF density reached a high level (ρ = 60/m or 72/m), SCNF became dominant, both at the border of the SRV and around the injection hole.TCRM was only sparsely distributed in HFN, but the number of these cracks was still far more than that of SCRM and TCNF.Figure23demonstrates the incremental accumulation of TCRM, SCRM, TCNF, and SCNF under different CPRs during fracturing process.In the rock formations with the same NF densities, the total number of cracks under S H : S h = 10:10 was less than that in the other two CPRs.The results suggested that as the CPR approached 1.0, it became more difficult for the fracturing treatment to form rocks.When ρ = 24/m, TCRM had the most value under all three CPRs.With an increase of NF density (ρ = 36/m), SCNF values grew gradually and reached roughly an equal level with TCRM value.When ρ ≥ 48/m, the SCNF value overtook TCRM value and became dominant among different crack types.The higher the NF density, the bigger of the gap between the value of SCNF and those of the other crack types.The value of TCNF also grew gradually with the increase in NF density.When ρ = 72/m, the SCNF value almost doubled TCRM value.In particular, NF density and confining pressure had little influence on the SCRM value.

Figure 22 .
Figure 22.Change of crack distribution with confining pressure ratio when NF density ranges from 24/m to 72/m and injection rate is 2.0 × 10 −6 m 3 /s.

Figure 22
Figure 22 demonstrates the influence of NF density on crack distribution under different CPRs.

Figure 24 .Table 6 .
Figure 24.Variation of total number of cracks with NF density under different confining pressure ratios and injection rate is 2.0 × 10 −6 m 3 /s.
Figure 24 reveals an exponential relationship between the total number of cracks and NF density.When NF density was relatively low ( ρ = 24/m or 36/m), different CPRs resulted in large differences in the total number of cracks, which reached the maximum level when SH: Sh = 10:5.The differences in the total number of cracks under different CPRs were gradually narrowed with an increase of NF density ( ρ ≥ 36/m).In other words, the variance of CPR had little influence on the change in the total number of cracks in rock formations with a high NF density.

Figure 24 .Table 6 .
Figure 24.Variation of total number of cracks with NF density under different confining pressure ratios and injection rate is 2.0 × 10 −6 m 3 /s.

Figure 24 .
Figure 24.Variation of total number of cracks with NF density under different confining pressure ratios and injection rate is 2.0 × 10 −6 m 3 /s.

Figure 25 .
Figure 25.Change of injection pressure with injection rate when NF density ranges from 24/m to 72/m and SH: Sh = 10:5.

Figure 25 .
Figure 25.Change of injection pressure with injection rate when NF density ranges from 24/m to 72/m and S H : S h = 10:5.

Figure 25 .
Figure 25.Change of injection pressure with injection rate when NF density ranges from 24/m to 72/m and SH: Sh = 10:5.

Figure 26 .Table 7 .
Figure 26.Variation of breakdown pressure with NF density under different injection rates and confining pressure ratio is S H : S h = 10:5.

Figure 27 .
Figure 27.Change of crack distribution with injection rate when NF density ranges from 24/m to 72/m and confining pressure ratio is S H : S h = 10:5.

Figure 28
Figure 28 demonstrates the incremental accumulation of TCRM, SCRM, TCNF, and SCNF under different injection rates.In rock formations with ρ = 24/m, TCRM always predominated, regardless of the changes in injection rate.When NF density increased from 24/m to 48/m or 60/m, larger injection rates led to more TCRM and SCNF, and SCNF dominated TCRM in value under relatively high injection rates.When NF density reached the maximum value (ρ = 72/m), SCNF almost doubled TCRM in value.Similarly, the changes of NF density and injection rate had little influence on the value of SCRM and TCNF.

Figure 27 .
Figure 27.Change of crack distribution with injection rate when NF density ranges from 24/m to 72/m and confining pressure ratio is SH: Sh = 10:5.

Figure 28 Figure 28 .
Figure 28 demonstrates the incremental accumulation of TCRM, SCRM, TCNF, and SCNF under different injection rates.In rock formations with ρ = 24/m, TCRM always predominated, regardless of the changes in injection rate.When NF density increased from 24/m to 48/m or 60/m, larger injection rates led to more TCRM and SCNF, and SCNF dominated TCRM in value under relatively high injection rates.When NF density reached the maximum value ( ρ = 72/m), SCNF almost doubled TCRM in value.Similarly, the changes of NF density and injection rate had little influence on the value of SCRM and TCNF.

Figure 29 .Table 8 . 6 0
Figure 29.Variation of the total number of cracks with NF density under different injection rates and confining pressure ratio is SH: Sh = 10:5.

Figure 29 .Table 8 .
Figure 29.Variation of the total number of cracks with NF density under different injection rates and confining pressure ratio is S H : S h = 10:5.

Figure 30 .
Figure 30.Change of injection pressure and crack distribution with particle size when NF density is 48/m and SH: Sh = 10:5.

Figure 30 .
Figure 30.Change of injection pressure and crack distribution with particle size when NF density is 48/m and S H : S h = 10:5.

Figure 31 .
Figure 31.Change in injection pressure and crack distribution with model size when NF density is 48/m and SH: Sh = 10:5.

Figure 31 .
Figure 31.Change in injection pressure and crack distribution with model size when NF density is 48/m and S H : S h = 10:5.

Table 2 .
Micro-properties used in smooth joint (SJ) contact for simulated specimens after calibration.

Table 3 .
Computational parameters of fluid properties.

Table 4 .
Fitting curve equations of the number of cracks for different NF densities.

Table 6 .
Fitting curve equations of total number of cracks under different confining pressure ratios.