Numerical Investigation of Influence of In-situ Stress Ratio, Injection Rate and Fluid Viscosity on Hydraulic Fracture Propagation Using a Distinct Element Approach

Numerical simulation is very useful for understanding the hydraulic fracturing mechanism. In this paper, we simulate the hydraulic fracturing using the distinct element approach, to investigate the effect of some critical parameters on hydraulic fracturing characteristics. The breakdown pressure obtained by the distinct element approach is consistent with the analytical solution. This indicates that the distinct element approach is feasible on modeling the hydraulic fracturing. We independently examine the influence of in-situ stress ratio, injection rate and fluid viscosity on hydraulic fracturing. We further emphasize the relationship between these three factors and their contributions to the hydraulic fracturing. With the increase of stress ratio, the fracture aperture increases almost linearly; with the increase of injection rate and fluid viscosity, the fracture aperture and breakdown pressure increase obviously. A low value of product of injection rate and fluid viscosity (i.e., Qµ) will lead to narrow fracture aperture, low breakdown pressure, and complex or dispersional hydraulic fractures. A high value of Qµ would lead wide fracture aperture, high breakdown pressure, and simple hydraulic fractures (e.g., straight or wing shape). With low viscosity fluid, the hydraulic fracture geometry is not sensitive to stress ratio, and thus becomes a complex fracture network.


Introduction
Hydraulic fracturing technology is widely used in the unconventional energy field [1][2][3].In recent years, hydraulic fracturing technology has played a critical role in the shale gas revolution [4,5].The in-situ stress ratio, injection rate, and fluid viscosity all have a crucial effect on the hydraulic fracturing results.King [6] reviewed the 30-year experience in hydraulic fracturing of shale in North America, and showed that the desired hydraulic-fracturing effect depends on not only rock mechanics in the shale reservoir, but also on in-situ stress state and operational parameters.Some researchers [7][8][9] have investigated the influence of injection rate on hydraulic fractures through laboratory experiments, which indicated that under the same experimental conditions, with the increase of injection rate, the branches of the hydraulic fracture would decrease, and the complex fracture network turns into a single straight fracture.Ishida et al. [10][11][12][13], Fan and Zhang [9], and Hou et al. [14], have performed laboratory experiments on the effect of fluid viscosity.Warpinski et al. [15] have analyzed the microseismic monitoring data in hydraulic fracturing sites.These studies indicated that under the condition of low fluid viscosity (such as the case of supercritical CO 2 ), the hydraulic fracture becomes flexed with multiple branches and develops a complex fracture network.With the increase of fluid viscosity (such as the case of a gel), the branches disappear and hydraulic fracture turns into a simple wing shape.We note that in these studies, only Ishida et al. [10][11][12][13] used the supercritical CO 2 and liquid CO 2 as fracturing fluid of which the viscosity is smaller than 1 mPa¨s.Supercritical CO 2 is important but difficult to use in hydraulic fracturing [16], thus, we need to perform systematic studies to enhance the theoretical basis for unconventional resource production.Meanwhile, the relationship between stress ratio, injection rate, and fluid viscosity should be further investigated and analyzed.
Numerical simulation [17][18][19][20][21][22][23][24][25][26][27][28] would avoid some of the randomness in laboratory experiments and could reflect the effect of critical parameters on hydraulic fracturing characteristics.Many works have been carried out to investigate the influence of some important parameters on hydraulic fracturing by numerical simulation.Wangen [29] established a two-dimensional numerical model based on the finite element method and further analyzed the bottom-hole pressure and the number of failed elements under different permeability conditions.It was found that under conditions of low permeability, the bottom-hole pressure rises more quickly and more elements fail.Carrier and Granet [30] developed a cohesive zone to govern the fracture propagation in the finite element analysis framework.Carrier and Granet's numerical results about injection pressure and fracture aperture are very consistent with analytical solutions under a zero fluid lag assumption.Zhang et al. [31] set up a 3D nonlinear fluid-mechanics coupling model, by using finite element software ABAQUS.The field data from the Daqing (in China) oilfield were used in numerical computations, and the evolution of the bottomhole pressure in the simulation was coincident with the corresponding data from the field measurements.Subsequently, the influence of some critical parameters on hydraulic fracturing characteristics were investigated.Using the finite difference software Flac 3D , Zhou and Hou [32] established a model in which the leak-off could be simulated in both rock matrix and fractures.They verify the modeling approach by comparing the borehole pressure between numerical and laboratory results.They further applied this approach on a tight gas sandstone reservoir from the Northern German Basin.In addition, some commercial software based on a continuum were launched for design and prediction of hydraulic fracturing, such as StimPlan, FRANC3D, and GOHFER.However, the relationships between stress ratio, injection rate, and fluid viscosity have not been comprehensively investigated, and the function of supercritical CO 2 with low fluid viscosity is still not clear.Meanwhile, previous studies conducted with discontinuum numerical method have rarely been reported.
In this paper, the Universal Distinct Element Code (UDEC) [33], a distinct element method (DEM) [34], was used to simulate hydraulic fracture initiation and propagation.The numerical solution obtained is consistent with the analytical solution of the breakdown pressure, which indicates that the distinct element method is feasible for modeling the hydraulic fracturing.Then, based on our numerical model, the influence of stress ratio, injection rate, and fluid viscosity on hydraulic fracturing are independently investgated by tuning one while retaining the other parameters.Finally, the relationships between the above three parameters are analyzed.

Algorithms Introduction
The DEM is a numerical simulation method which is based on a discontinuum approach.It was put forward by Cundall [34] in 1971, and now has become an effective method for dealing with geotechnical problems.The DEM allows finite displacements and rotations of discrete bodies including complete detachment, and recognizes new contacts automatically as the calculation progresses [35].With the DEM software Particle Flow Code (PFC), some works have been carried out to simulate the hydraulic fracturing characteristics [36][37][38][39].However, there are some disadvantages in simulating hydraulic fracturing using PFC.The PFC assumes that each element acts as a rounded particle; thus, the PFC could not well describe the contact relationships between particles.In other words, the PFC does not satisfy the failure criterion, for example, the Mohr-Coulumb failure criterion.In addition, the PFC would generate huge amounts of particles when simulating hydraulic fracturing on a large scale, which results in an extremely huge computational cost.In contrast, the UDEC is immune to these disadvantages because it assumes that each particle is polygonal.The particle, or block, in the UDEC is impermeable, and the fluid flows only in the contacts between blocks, which is consistent with the physical characteristics of the extremely low permeability of the shale matrix.Using UDEC, Choi [40] proposed several methods for obtaining shut-in pressure from the ambiguous shut-in curves, under various remote stress regimes and various rock properties.Zangeneh et al. [41] investigated the interaction between the natural fractures and hydraulic fracture, and the results show that the natural fractures influence significantly the size, orientation, and path of the hydraulic fracture.
In the UDEC, the research object is discretized into rigid or deformable blocks by one group or several groups of discontinuity.If the blocks are deformable, they will be subdivided into a mesh of finite-difference elements, and each element responds according to a prescribed linear or nonlinear stress-strain law.These rigid or deformable blocks are considered to be impermeable to the flow of fluid restricted to the contacts.For edge-to-edge contacts in fluid-mechanical analyses, the calculation of flow rate q obeys the cubic law [33]: where q is the flow volume between contacts per unit time, a is the joint hydraulic aperture, µ is the dynamic viscosity of fluid, ∆p is the pressure difference between the domains, and L is the contact length.
The contact conductivity depends on deformation and displacement, which is in turn influenced by the fluid pressure.Figure 1 shows the relationship between aperture and normal stress.Given the hydraulic aperture at zero normal stress a zero , the minimum allowable joint hydraulic aperture (residual aperture) a res , and the maximum allowable joint hydraulic aperture a max, we can simulate the hydraulic fracture.When the hydraulic aperture exceeds a zero , or, the effective stress acting on the joint is 0, the hydraulic fracture is assumed to be open and will propagate along the contacts between blocks.

Model Setup and Parameters Acquisition
A two dimensional model was built with UDEC, to investigate the hydraulic fracturing characteristics of laboratory-scale shale.According to laboratory hydraulic fracturing experiments [7,14,42,43], the model dimension was confirmed to be 300 mm ˆ300 mm.There is a circular hole that is 20 mm in diameter at the center of the model (Figure 2), from which the fluid was injected at a constant rate.The fluid is considered as incompressible with transient flow.According to references [40,41], an irregular distribution of elements leads to an irregular hydraulic fracture.This would disturb the influence of critical factors we investigated here.To eliminate the influence of the irregular shape, the model is discretized into 900 squares with 10 mm side length.In view of the elastic-brittle mechanical characteristic of rock, the blocks are considered to be elastic to reflect the deformation characteristic of the model.The contacts are considered to be elastic-plastic with Mohr-Coulomb failure criterion.When interacting with natural fracture, hydraulic fracture will turn to be asymmetrical and complicated.Therefore, the influence of natural fractures on hydraulic fracturing was not considered in our model.As the research object is intact rock, the contact strength parameters are assumed to be high (say close to the intact rock strength), which makes the mechanical response of the numerical model close to the intact shale.Hydraulic fracture will propagate along the contacts between blocks, and the discrete block may deform and move.To investigate the laboratory mechanical and physical properties and provide accurate parameters for numerical modeling of hydraulic fracturing, shale specimens were collected from the outcrops of the Longmaxi Formation in Shizhu County (in China).The values of some parameters in Table 1 are rounded based on above laboratory experiments.The contact strength parameters are set to be high to make the mechanical response of the numerical model close to the intact rock.The normal stiffness of contacts, k n , was calculated [44]: where k is bulk modulus, g is shear modulus, and ∆Z min is the smallest width of the zone adjoining the contact.
The operational parameters, injection rate, and fluid viscosity are shown in Table 2.The stress boundary condition in the numerical simulation is that the major horizontal stress was set equal to 15 MPa (σ H = 15 MPa), and the minor horizontal stress was set equal to 10 MPa (σ h = 10 MPa).

Simulation Schemes
Based on the above configurations of the model, the initiation and propagation of hydraulic fracture for the numerical model was analyzed using UDEC, and the obtained numerical solution of the breakdown pressure was compared with the analytical solution.Then, the independent influence of stress ratio, injection rate, and fluid viscosity on hydraulic fracturing characteristics was investigated within reasonable ranges.Finally, the relation between above three parameters was analyzed.

Initiation and Propagation of Hydraulic Fracture
When fluid is first injected in the model, injection pressure rises linearly with time, as shown in Figure 3.When the running time equaled 1.5 s and the injection pressure reached 9.81 MPa, two hydraulic fractures initiated simultaneously on both sides of the hole (the grids were hidden).This indicates that the breakdown pressure should be 9.81 MPa.With the failure of contacts between blocks, the injection pressure was suddenly released and fluid leakage occurred because the aperture grew, which leads to an obvious pressure drop at 2 s.At 6 s (Figure 4), the hydraulic fracture further propagated.Every extension of the hydraulic fracture coupled with a corresponding pressure drop.The breakdown pressure P b was calculated [45]: where α is the poroelastic coefficient, υ is the Poisson ratio, and T is the rock tensile strength.We test α = 0.15 and 0.24 to calculate the breakdown pressure, respectively.These two cases (i.e., α = 0.15 and 0.24) can be considered as analytical solution, which were compared with our numerical solution to verify the reliability of the numerical simulation (Figure 5).Stress ratio K was defined to be σ H /σ h , where σ H is the major horizontal stress and σ h is the minor horizontal stress.The relation between breakdown pressure and K was showed in Figure 5, where we kept σ h = 10 MPa for simplicity.As shown in Figure 5, with the increase of stress ratio, the analytical and numerical solutions of breakdown pressure decreased, and the numerical solution decreases much more slowly.It is worth noting that when K is between 1.3 and 1.6, the numerical solution is close to the analytical solution.Therefore, the stress ratio for the basic model was set to 1.5, to assure the accuracy and reliability of the numerical simulation.

Influence of Stress Ratio
In this section, the influence of stress ratio on hydraulic fracturing characteristics was investigated on the condition that the running time equaled to 15 s and other parameters in the basic model remained unchanged.As shown in Figure 6, when the stress ratio K = 1.1~1.2(σ h is kept 10 MPa), the fractures were divergent and there were small branch cracks on both sides of the main hydraulic fracture.When the stress ratio K = 1.3, the main hydraulic fracture propagated, and most branch cracks disappeared.When the stress ratio K ě 1.4,only the main hydraulic fracture existed, and the branch cracks totally disappeared.The main reason for this result is that hydraulic fracture usually propagates along the direction perpendicular to minimum principal stress.When the stress ratio K is close to 1, i.e., the stress field is nearly isotropic, the branches of hydraulic fracture will appear.Whereas, with the increase of stress ratio K, the stress anisotropy is more and more significant; thus, those branch cracks would disappear.Figure 7 shows the relation between hydraulic fracture aperture and stress ratio.Obviously, the aperture increases with increasing K.The increase of stress ratio would cause the decrease of breakdown pressure; thus, the earlier generated hydraulic fracture would lead to a much wider aperture for the same amount of time.

Influence of Injection Rate
In this section, the influence of injection rate on hydraulic fracturing characteristics was investigated on the condition that running time equaled to 15 s and the other parameters in the basic model remained unchanged.As shown in Figure 8, when Q = 8 ˆ10 ´6 m 3 /s/m, the branch crack initiated.When Q ě 9 ˆ10 ´6 m 3 /s/m, branch cracks disappear and the main hydraulic fracture extends significantly with an increasing injection rate.The lower injection rate results in the more lateral fluid leak-off, which makes a much smaller branch crack.As shown in Figure 9a, the hydraulic fracture aperture is also affected by injection rate, and the aperture would increase obviously with increasing injection rate.Figure 9b shows that with the increasing injection rate, the breakdown pressure increases almost linearly.When the injection rate is high, the loading condition would approach to dynamic procedure; thus, the breakdown pressure needed for hydraulic fracture initiation is higher.This phenomenon is consistent with the rock uniaxial compression tests.
Some researchers have investigated the influence of injection rate on hydraulic fractures taking water [7], viscasil [8], and gel [9] as fracturing fluid.These studies indicated similar conclusions that under the same experiment conditions, with the increase of injection rate, the branches of the hydraulic fracture decrease, and the complex fracture network turns into a single straight fracture.

Influence of Fluid Viscosity
In this section, the influence of fluid viscosity on hydraulic fracturing characteristics was investigated on the condition that running time equaled to 10 s and other parameters in the basic model remained unchanged.As shown in Figures 10 and 11a, the length and aperture of hydraulic fracture increases significantly with increasing viscosity.The main reason is that the higher fluid viscosity would diminish the lateral fluid leak-off, which noticeably elevates the fluid volume needed for generating straight fracture.With the increasing fluid viscosity, breakdown pressure increases almost linearly, as shown in Figure 11b.This is because higher viscosity diminishes the fluid leak-off, which is equivalent to elevating the fluid volume needed for fracturing within the same time.According to the analysis of Section 3.3, the above reason leads to higher breakdown pressure.Some researchers, through laboratory hydraulic fracturing experiments [8][9][10][11][12][13][14] and site microseismic monitoring data [15], arrived at the similar conclusions that under the condition of low fluid viscosity (such as the case of supercritical CO 2 ), the hydraulic fracture becomes flexed and with multiple branches, and develops a complex fracture network; with the increase of fluid viscosity (such as the case of a gel), the branches disappear and the hydraulic fracture becomes a simple wing shape.
According to [12,13], the viscosity of supercritical CO 2 and liquid CO 2 are 0.05 mPa¨s and 0.1 mPa¨s, respectively.Under the same condition, we simulated the effect of these two fluids on hydraulic fracturing, and no fracture appeared.This result indicates that when the injection rate is small enough, a low viscosity fluid could not break the rock.Based on the above numerical simulation and previous laboratory experiment results [9,14], as shown in Figure 12, we could know that when the fracturing fluid viscosity is too low, the injection rate could be elevated to produce hydraulic fracture.We verified the reliability of the numerical model, through comparison between numerical simulation, analytical solution, and experiment results shown in Sections 3.1-3.4.
It is obvious that when the flow rate q or fluid viscosity µ is increasing, the influence is equal on the right side of Equation ( 4), that is, either would lead to an increase of both the aperture a and the pressure ∆p.Because there is a positive correlation between injection rate Q (a parameter in software UDEC) and flow rate q, we could replace q with Q.The product Qµ is considered to be a whole in this section, to investigate its influence on hydraulic fracturing characteristic.
In Section 3.4, we did see the hydraulic fracture using the fluid viscosity µ = 0.05 mPa¨s and 0.1 mPa¨s.Thus, in this section, we take Q = 100 ˆ10 ´6 and 50 ˆ10 ´6 m 3 /s/m to have Qµ = 5 ˆ10 ´6 mN; meanwhile, when fixing the fluid viscosity at 1 and 2 mPa¨s, we adjust the injection rate to have Qµ = 10 ˆ10 ´6 mN, as shown in Table 3.For the convenience of comparison, we retain the stress ratio at K = 1.5 and the time at 10 s to verify the influence of the Qµ on the hydraulic fracturing.Figures 13 and 14 show that the hydraulic fracture becomes flexed and complex with multiple branches when Qµ is small, and the corresponding fracture aperture as well as the breakdown pressure are small.In contrast, the branches disappear and hydraulic fracture turns to be a simple wing shape.Therefore, a small value of Qµ would be helpful for generating complex fracturing network.

Relationship between Stress Ratio and Fluid Viscosity
As shown in Section 4.1, we can obtain hydraulic fractures if we increase the injection rate using supercritical CO 2 with fluid viscosity of 0.05 mPa¨s.Thus, we verify the sensitivity of the low viscosity fluid on the stress ratio by only tuning the stress ratio while retaining the other parameters.
With the stress ratio of K = 1.1, the hydraulic fracturing does not occur; thus, the results are not shown in Figures 15 and 16    As shown in Figure 15, for low viscosity fluid (e.g., supercritical CO 2 ), the geometry of the hydraulic fractures are not sensitive to stress ratio, and still turns to be flexed and complex, which is quite different to those shown in Figure 6.This result indicates that in the formation with remarkable anisotropy stress, taking supercritical CO 2 as fracturing fluid is in favor of hydraulic fracturing.The stress ratio still influences the breakdown pressure of the rock, as shown in Figure 16b).

Conclusions
Based on the shale mechanical parameters from laboratory tests, we set up the basic configuration of a numerical model, which was used to verify the reliability of the distinct element code UDEC on simulating the hydraulic fracturing.We mainly discussed the influence of critical parameters on hydraulic fracturing.The results are summarized as follows: (1) The increasing injection rate and fluid viscosity have the same influence on hydraulic fracturing characteristics.With the increasing injection rate and fluid viscosity, fracture aperture and breakdown pressure obviously increase.(2) The hydraulic fracture becomes flexed and complex with multiple branches when the value of Qµ is low, and the corresponding fracture aperture as well as the breakdown pressure are small.In contrast, the branches disappear and hydraulic fracture becomes a simple wing shape, and the corresponding fracture aperture as well as the breakdown pressure are high.(3) With low viscosity fluid, the hydraulic fracture geometry is not sensitive to stress ratio, and becomes a complex fracture network.

Figure 1 .
Figure 1.The relationship between fracture aperture and normal stress.

Figure 3 .
Figure 3.The hydraulic fracturing characteristics at the time of 2 s.(a) fracture geometry; (b) injection pressure curve; (c) distribution of pore pressure.

Figure 5 .
Figure 5.Comparison of breakdown pressure under different stress ratios between the analytical and numerical solutions.

Figure 7 .
Figure 7.The relationship between maximum aperture and stress ratio.

9 .
The relationship between hydraulic fracturing characteristics and injection rate.(a) The relationship between maximum aperture and injection rate; (b) The relationship between breakdown pressure and injection rate.

Figure 11 .
Figure 11.The between hydraulic fracturing characteristics and injection fluid viscosity.(a) The relationship between maximum aperture and fluid viscosity; (b) The relationship between breakdown pressure and fluid viscosity.

Figure 12 .
Figure 12.The relationship between injection rate and fluid viscosity.

4 .
Relationship between in-situ Stress Ratio, Injection Rate and Fluid Viscosity4.1.Relationship between Injection Rate and Fluid ViscosityRearranging Equation (1), we obtain:

Figure 13 .
Figure 13.The relationship between hydraulic fracture geometry and Qµ.

Figure 14 .
Figure 14.The between Qµ and hydraulic fracturing characteristics.(a) relationship between Qµ and the fracture aperture; (b) relationship between Qµ and the breakdown pressure. .

Figure 16 .
Figure 16.The relationship between ratio and hydraulic fracturing characteristics.(a) relationship between stress ratio and the fracture aperture; (b) relationship between stress ratio and the breakdown pressure.

Table 1 .
Block and contact parameters in numerical model.

Table 3 .
Different combination of the injection rate and viscosity.