Numerical Analysis on the Optimization of Hydraulic Fracture Networks

The clear understanding of hydraulic fracture network complexity and the optimization of fracture network configuration are important to the hydraulic fracturing treatment of shale gas reservoirs. For the prediction of hydraulic fracture network configuration, one of the problems is the accurate representation of natural fractures. In this work, a real natural fracture network is reconstructed from shale samples. Moreover, a virtual fracture system is proposed to simulate the large number of small fractures that are difficult to identify. A numerical model based on the displacement discontinuity method is developed to simulate the fluid-rock coupling system. A dimensionless stress difference that is normalized by rock strength is proposed to quantify the anisotropy of crustal stress. The hydraulic fracturing processes under different stress conditions are simulated. The most complex fracture configurations are obtained when the maximum principle stress direction is perpendicular to the principle natural fracture direction. In contrast, the worst results are obtained when these two directions are parallel to each other. Moreover, the side effects of the unfavorable geological conditions caused by crustal stress anisotropy can be partly suppressed by increasing the viscous effect of the fluid.


Introduction
The dense spacing hydraulic fracture network is important for shale gas reservoirs [1,2].Generally, there are two kinds of factors related to the hydraulic fracturing process.One kind is the geological parameters, such as the rock properties and crustal stress conditions.Another kind is the operating parameters, such as the injection rate and fluid viscosity.Given the geological conditions, it is important to understand how fractures propagate and to find the optimal operating parameters under these geological conditions.
Fracturing treatment designs rely heavily on imprecise estimates of the stimulated reservoir volume derived from micro-seismic observations and through a highly inefficient trial-and-error approach [3].However, the details of fracture networks cannot be measured by these observation methods.The fracture propagation has also been observed in laboratory [4][5][6].However, it is difficult to observe fracture network propagation in the laboratory and to know how laboratory experiments relate to the reservoir scale [7].Compared with the above methods, numerical modelling is much more flexible.Numerical models can be established according to field conditions and all the details can be exported for analysis.Therefore, numerical modeling is an important tool for engineers to predict the geometry of fracture networks [8].
Numerical models should capture all the essential elements so that a simulation reasonably represents the real process [7].However, the rigorous simulation of all the aspects is very challenging.Many early stage models, such as the PK model [9], PKN model [10], KGD model [11], pseudo-3D models and planar-3D models [12,13], can only be used in simulating the propagation of fractures with highly idealized geometries.Models for the simulation of complex fracture networks are still in development in recent years.These models can roughly be classified into three kinds.The first kind is the models which need grid meshing in the whole region, such as the finite difference method [14], finite element method [15], extended finite element method [16][17][18] and discrete element method [19,20], etc.There are many softwares based on these methods, such as UDEC/3DEC, PFC, RFPA, FLAC, etc.However, while these methods are easy to use, they are far from perfect.As the simulation domain is often several orders bigger than the aperture of fractures, a deformation that is considered as small "noise" in a solid solver may induce dramatic oscillations of fluid pressure in a flow solver [19].As a result, these methods are ill-conditioned for the fluid-rock coupling system of hydraulic fracturing.Another problem is the complexity of grid re-meshing, which is necessary for the accurate simulation of fracture trajectory [20].The second kind is meshless methods.Different from traditional methods, meshless methods [21][22][23] are free from meshing and have the capability to deal with moving boundary conditions.Meshless methods are flexible and computationally efficient in simulating fracture propagation, but have not been in used in simulating the complex fracture network propagation in hydraulic fracturing.
The third kind of fracture network simulation model is the models based on Displacement Discontinuity Method (DDM) [24][25][26][27], which is a kind of boundary element method.In DDM-based models, grid meshing is needed only for fractures.The propagation of a fracture network can be simulated within less CPU time because the grid number is far less than in the above models.Fracture aperture can be calculated with higher precision because the stress interaction between fractures is directly calculated by analytical solutions.Moreover, when fractures extend, the grid can be adjusted by simply adding new fracture elements.Based on DDM, Zhang and Jeffrey [28] investigated the injection type on the creation of complex fracture network and it was found that the equal rate injection is better than the equal pressure injection.McClure and Horne [29] simulated the propagation of complex fracture network and showed that fractures may propagate through a mixture of the opening and sliding of preexisting fractures and the propagation of new fractures.Zhang et al. [30] investigated the stability of hydraulic fracture network propagation and found that the localization instability can be suppressed by fluid viscosity.To sum up, the DDM-based models are appropriate for the simulation of complex fracture networks and have been used by more and more researchers [3,7,8,[28][29][30][31][32][33].
The Natural Fracture (NF) characterization is also a critical foundation.NFs often lead to the formation of complex fracture geometries by altering the propagating direction of Hydraulic Fracture (HF).However, it's difficult to obtain the geometry of NFs.For theoretical purposes, idealized NF geometries are acceptable [34], but for engineering purposes, complex fracture networks are required.The Discrete Fracture Network (DFN) model is always used to reconstruct complex NF networks.With the distributions of fracture length, orientation, aperture, etc., and the fracture density, random NF networks can be created [20].The DFN model has been used by many researchers [3,7,8,20,32].For example, Riahi and Damjanac [20] investigated the effects of NF connectivity, injection pressure, etc., on the complexity of hydraulic fracture network.Kresse et al. [32] investigated the propagation of fracture networks when complex natural fractures are created.McClure and Horne [7] investigated how shear stimulation can affect the process of fluid propagation.In summary, it's a useful method to reconstruct the complexity of NF network by random distributions.However, no matter how many distributions are used, the randomly created NF networks are never totally equivalent to the real fracture network.
In this work, the NF network is reconstructed from shale samples.Considering the characteristics of the real NF network, a virtual fracture system is proposed to deal with the small fractures that are difficult to identify from shale samples.A DDM-based numerical model is developed to investigate how hydraulic fracture networks propagate under different stress conditions and how to get complex hydraulic fracture network by adjusting the viscous effect of the fluid.

Basic Equations
The fluid-rock coupling process of hydraulic fracturing is simulated based on DDM.The following assumptions are used: the domain of the rock matrix is infinite and the rock matrix is homogeneous, isotropous and linearly elastic [35].The rock matrix is impermeable.The fluid injected is Newtonian, single phase and laminar [7,28,36].
As a kind of boundary element method, the theoretical solution of elastic mechanics [35] is directly used in DDM.Given a fracture system with N fractures, the stress induced by the opening and sliding of these fractures is calculated by [28]: where x = (x, y) is the coordinate, w is the normal displacement discontinuity, which equals to the fracture aperture, v is the shear displacement discontinuity, which equals to the sliding distance of fracture surface, lr is the length of fracture r, Gij are the hyper singular Green's functions, which are proportional to the plane strain Young's modulus [28], σn is the normal stress and τs is the shear stress, obeying Coulomb's frictional law characterized by the friction coefficient λ, which limits the shear stress by: that can act in parts of fractures that are in contact, but vanish along the separated parts.Along the opened fracture portions with fluid pressure pf, we have: K is the three dimensional correction coefficient.Using the parameters given by Wu and Olson [37], the correction coefficient K proposed by Olson [38] can be written as: where h is the limite layer thickness perpendicular to the simulation plain, d is the distance between points x 1 and x 2 .
Volumetric flux q of fluid in fracture is described by: where μ is fluid dynamic viscosity, w0 , which is a reflection of the fracture surface roughness, is the initial hydraulic aperture, w is the aperture, The global mass balance requires: where vf is fluid volume, L is fracture segment length, q is the volumetric flux of fluid.The fracture growth is based on the maximum hoop stress criterion, with the maximum mixed-mode intensity factor reaching a critical value: where KI and KII are stress intensity factors, KIC is tensile mode fracture toughness, θ is the fracture propagation direction relative to the current fracture orientation and satisfies: The fractures are meshed with constant displacement elements.The hydraulic fracture problem is solved simultaneously including the effects of viscous fluid flow and coupled rock deformation.

Natural Fracture Network
NF networks are measured from shale samples, which are sampled from the outcrops of the Longmaxi Formation in South China.As the fracture aperture is always no bigger than 100 μm, it's impossible to directly obtain the fracture trajectory.Therefore, water was sprinkled over the shale sample.Some water enters NFs so that the rock surface is darker near NFs.The surfaces of the shale samples are shown in Figure 1a-c Two characteristics are important to our modelling: (1) The complete identifying of NF is impossible.By contrast, many small fractures are not clearly shown.(2) There are two orthogonal sets of fractures.The length is longer and the connectivity is better for the horizontal set.

Virtual Fracture System
There are large number of small fractures that are difficult to be identified.As shown in Figure 2, there are many small fractures that are linked to the main fractures.In fact, there are many other fractures that cannot be seen in the photos.Although the aperture of these fractures is small and the length is short, the strength of these fractures is smaller than the rock matrix.Therefore, the effect of these small fractures is not ignorable.
A virtual fracture system is introduced to simulate these small fractures.First, when a HF encounters a NF, a virtual fracture is added at the other side of the NF.Second, virtual fractures are added randomly to some fracture elements.The methodologies are illustrated in Figure 3.During the simulation, the virtual fractures are checked under the in situ stress conditions.If propagation can occur for a virtual fracture, it means that the in situ stress is so big that the rock may break at the location of the virtual fracture.Under this circumstance, a new fracture element would be added.Otherwise, the small fractures can be treated as if they do not exist and the rock matrix can be treated as homogeneous.This is what the term "virtual" refers to.
Virtual fractures are useful in simulating the creation of complex fracture networks.When virtual fractures extend, new fracture branches are formed and thus the fracture network becomes more complex.The interaction between HF and NF can also be simulated conveniently using the virtual fracture system.The HF will across the NF when the virtual fracture element at the crossing point can extend.Moreover, as no stress is induced by the virtual fractures, no extra computationally demanding is needed.The using of virtual fracture in the modelling will be introduced in Section 2.4.

The Work Flow
The most important process of the modelling is the calculating of pf, w, v and vf, where pf is fluid pressure, w is aperture, v is sliding distance, vf is fluid volume.A Gauss-Seidel iteration method is used to calculate these values.The balance values of each element are calculated by finding the values which satisfy the Equations ( 2), ( 3), ( 5) and ( 6).The most time-consuming process is the calculating of in situ stress by Equation (1), which is required for solving these balance values and must be updated in each iteration.As the in situ stress is related to all fractures, the computational demands increase rapidly with the element number.Generally, more than 95% of CPU time is consumed in this process, which is shown as the red dashed box in Figure 4.It must be noted that the virtual fractures are not incorporated in this process.Therefore, although there are large number of virtual fractures, the effectiveness of the solution is not affected.
After pf, w, v and vf are updated, the virtual fractures and the ordinary fractures are checked individually to see if fractures could propagate.Specially, for the virtual fractures, as shown in the blue dashed box in Figure 4, we calculate the in situ stress and then the stress intensity factor.If the stress intensity factor is bigger than the fracture toughness, a new ordinary fracture element is added and the corresponding virtual fracture is removed.

Against Theoretical Solution
The solving of fracture aperture is a most important foundation.Given the fluid pressure, the profile of the aperture w of a single fracture in infinite rock matrix follows [35]: where a is the fracture half length, pf is the net fluid pressure, x is the distance to the center of the fracture, G is the shear modulus, νPo is the Poisson's ratio.The numerical model is validated against the analytical solution.A single fracture is simulated with the parameters listed in Table 1.The aperture profile of the fracture is compared with the analytical solution as shown in Figure 5.The numerical results agree well with the analytical one, even when there are only ten grids.
The precisely solving of aperture is especially important for the tip element because the intensity factor is calculated by the displacements of the tip element [39]: where κ 3 4 , w and v are the normal and shear displacements of the tip element respectively, l is the length of the tip element.Table 2 shows the precision of the displacements of tip element under different grid numbers.The error of the numerical simulation is significant.The reason is the using of the constant displacement element.Fortunately, the error is almost a constant.As shown in Table 2, the computed aperture is approximately 1.25 times bigger than the analytical value.Therefore, a correction coefficient is introduced, i.e.: where C = 1/1.25 is the correction coefficient.

Against Numerical Modelling
We refer for comparison to the model of Wu and Olson [8] for the fracture configurations of a single horizontal wellbore with two initial fractures (Figure 6a) and two horizontal wellbores with one fractures each (Figure 6b).These are the most classic configurations and have been simulated in many works [3,8,32,37,40,41].The correct simulation of these configurations requires the precisely solving of the stress field and fluid pressure, and the correct calculating of the stress intensity factors of fracture tips.These two configurations are simulated with the parameters listed in Table 3, which is the same to that of Wu and Olson [8].The results are shown in Figure 7.For the initial configuration as shown in Figure 6a, the two fractures propagate away from each other while for the initial configuration as shown in Figure 6b, the two fractures grow toward each other.In summary, our model compares favorably with previous models.

Numerical Setting
As shown in Figure 8, fluid is injected from the center of the region with constant flow rate qinj.As fluid must be injected in fractures, two small fractures as the red "+" shown in Figure 8 are added.Therefore, the initial fractures may propagate simultaneously along both x and y directions.This is helpful in reducing the influence of initial orientation on final configuration.
The rock matrix is assumed to be infinite.Therefore, no boundary conditions are needed.The maximum principle stress direction of far field stress is represented by ψ.The default input parameters are listed in Table 4, which is used in the following sections if there is no special statement.The simulation is ended when any HF reaches the boundary, i.e., the red dashed box in Figure 8.

Quantification of Stress Anisotropy and Fluid Viscosity
Generally, the stress anisotropy can be quantified by both the stress ratio and stress difference.In this part, we will first check which is more closely related to the propagation of fracture network.The effect of stress ratio is first checked.The fracturing processes under same stress ratio but different stress amplitudes are simulated.The results are shown in Figure 9, which shows different fracture configurations under same stress ratio.By contrast, for the fracturing processes under same stress difference, the final fracture configurations are almost equal as shown in Figure 10.These results indicate that stress difference is better in quantifying the stress anisotropy that related to fracture network configuration.However, it's not a good idea to directly use the stress difference to quantify the stress anisotropy because the fracture propagation are related not only to the stress condition, but also to the rock properties.The influence of rock properties can be partly removed by using the dimensionless stress difference which is normalized by rock strength.Considering the initiation fluid pressure can be used as a characterization of rock strength, the stress difference will be normalized by in this paper.Here refers to the fluid pressure above which hydraulic fracture begins to propagate.As discussed by Zhang et al. [30], is proportional to √ ⁄ , where s is the fracture spacing, then we define: where SA is the dimensionless stress difference, σ and σ are the two principle stresses, s is the average fracturing spacing of NF network.Specifically, when the principle stress directions are parallel to coordinate axis, the dimensionless stress difference can be written as: where σ and σ are the principle stresses along x and y directions respectively.
The viscous effect of fluid is an important operating parameter.Normalized by rock strength, a dimensionless number M is proposed to quantify the viscous effect of fluid by Zhang et al. [30]: where q is flow rate, μ is fluid dynamic viscosity, E is the Young's modulus.
Although fracture network propagation is related to many factors, including ∆σ, , q, μ and E, the effects of these factors are not independent.By contrast, the effects of these factors can be represented by two dimensionless numbers, i.e., dimensionless stress anisotropy SA and dimensionless viscous effect M as defined in the Equations ( 14) and ( 16).
To check the effectiveness of SA and M, we simulated the fracturing processes under same value of SA and M but different values of ∆σ, and μ as listed in Table 5.The simulation results are shown in Figure 11.Although all these parameters are different, the fracture configurations are almost equal as the dimensionless numbers are equal.These results indicate that stress anisotropy and fluid viscosity can be clearly quantified by SA and M, respectively.

Effect of Stress Anisotropy
First, we simulated the fracturing processes under different stress differences.All the other parameters are equal and are listed in Table 4.The results are shown in Figure 12.It's clear that stress difference has significant effect on fracture configuration.When 3.3, i.e., the maximum principle stress direction is parallel to x axis, HF propagates mainly along x axis.With the increase of σ , the propagation distance along y axis increases.Thus, the fracture network is more complex and the affected area is bigger when σ is bigger.The most complex fracture network is obtained when 0.66, under which circumstance the propagation is dominated by two opposite mechanisms.First, fractures tend to propagates along the principle NF direction, which is roughly parallel to x axis.Second, fractures tend to propagate along the maximum principle stress direction, which is parallel to y axis.The optimal results are obtained when the two opposite mechanisms are balanced by each other.Second, the fracturing processes under different values of ψ are simulated, where ψ is the angle between maximum principle stress direction and the positive direction of x axis as shown in Figure 8. ∆σ σ σ 9 MPa is used in the simulation.All the other parameters for these simulations are equal and are listed in Table 4.The simulation results are shown in Figure 13.The fracture configurations vary significantly with ψ.The propagation distance is further along the maximum principle stress direction.This is favorable for the creation of complex fracture network when the maximum principle stress direction is orthogonal to the principle NF direction.
The affected area is exported to investigate how stress anisotropy affects the fracturing result.The affected area is defined as the area of the region which has experienced a fluid pressure incensement due to the fluid injection [20].Specifically, a fracture is treated as affected when the fluid pressure incensement is bigger than 0.1 MPa in this paper.The variation of affected area with ψ and SA are shown in Figure 14.The affected area decreases with the increasing of SA and reaches its maximum value when the maximum principle stress direction is orthogonal to the principle NF direction.The variation of the average sliding distance of fracture in affected region with ψ and SA is shown in Figure 15.As the sliding of a fracture is caused by the shear stress along the fracture, the sliding distance vary significantly with the crustal stress direction.The sliding distance reaches its maximum value when ψ 45°, when the shear stress reaches the maximum value.However, the variation of sliding distance with ψ is significant only when the value of SA is big.Otherwise, if SA is very small, the sliding is also caused by the induced stress of other fractures, which never vary with the stress direction.

Effect of Fluid Viscosity
From the analysis of last section, we know the stress anisotropy has significant effect on the fracture network configuration.The problem is, it's hard to alter the stress condition during fracturing operation.Therefore, when the stress condition is unfavorable for fracturing treatment, the method to reduce the influence of stress anisotropy is very important.As the viscous effect of fluid is one of the most important operating parameters during fracturing treatment, the variation of affected area with SA and M is investigated.
The fracture network propagation under different values of SA and M are simulated.The value of SA is set by adjusting the value of ∆σ σ σ , the value of M is set by adjusting μ.All the other parameters are listed in Table 4.The affected area is plotted against SA and M in Figure 16, which shows a very clear tendency.When the injection rate is very low or viscous coefficient is small, i.e., the value of M is small, the fracturing result is good only when 0.6.By contrast, the optimal fracturing result can be obtained in a wider range of SAxy when the value of M is bigger.As the dimensionless number M is proportional to both q and μ, these results indicate that, if the crustal stress condition is unfavorable for fracturing treatment, the influence of stress anisotropy can be partly suppressed by increasing the injection rate or fluid viscosity.

Conclusions
A natural fracture network can be well represented by identifying the fracture trajectory directly from rock samples along with a virtual fracture system.The virtual fracture system is necessary because it's impossible to get all the details of NFs.The propagation of virtual fracture can be conveniently simulated within the framework of DDM without the loss of computational efficiency.The virtual fracture system is useful in simulating the HF-NF interactions and is instructive to other numerical methods.
The stress anisotropy can be quantified by defining a dimensionless stress difference SA, which is normalized by the strength of rock matrix.Given the certain value of SA and M, the fracture network has a certain configuration, which is not related to the injection rate, fluid viscosity, stress difference, rock strength, etc.
The propagation direction is dominated mainly by two mechanisms: the NFs and stress condition.As fractures tend to propagate along NFs and the direction of maximum principle stress, the complex fracture network can be formed when these two directions are orthogonal to each other and are balanced by each other.Under this circumstance, the leading propagation direction of HF is not exist.The key point of creating complex fracture network is the eliminating of the leading direction of HF propagation.
The side effect of stress anisotropy on the formation of complex fracture network can be weakened by increasing the viscous effect of fluid.Therefore, if the stress condition is unfavorable to fracturing operation, it is better to increase the injection rate or fluid viscosity.

Figure 1 .
Figure 1.The natural fracture networks created from the rock sample: (a-c) are the surfaces of the rock sample that wetted by water; (d-f) are the main fractures obtained from the rock sample.

Figure 2 .
Figure 2. The small fractures around main fractures.

Figure 3 .
Figure 3.The setting of virtual fractures.The blue line represents the NF, the black line represents the HF, and the small red lines represent the virtual fracture elements.

Figure 4 .
Figure 4.The flowchart of the model.

Figure 5 .
Figure 5.The comparison of numerical modelling with analytical solution.

Figure 6 .Figure 7 .
Figure 6.Two initiation fractures in (a) a horizontal wellbore and (b) two horizontal wellbores.Blue lines represent horizontal wellbores, and red lines represent initial fractures.

Figure 14 .
Figure 14.The variation of affected area with ψ and SA.The result of each simulation is represented by a black point.The smoothed value is represented by the colored surface.

Figure 15 .
Figure 15.The variation of shear displacements with ψ and SA.The result of each simulation is represented by a black point.The smoothed value is represented by the colored surface.

Figure 16 .
Figure 16.The variation of affected area with SA and M. The result of each simulation is represented by a black point.The affected area smoothed from the data points is represented by the colored surface.

Table 1 .
Input parameters for the calculation of fracture aperture profile.

Table 2 .
The comparison of numerical modelling with analytical solution for the aperture of tip element.Here x is the coordinate of the center of tip element.

Table 4 .
The default values of the input parameters.