Numerical Study on the Formation of Shear Fracture Network

Shear fracture network is important to the hydraulic fracturing treatment of a shale gas reservoir. In this paper, the formation of shear fracture network is investigated by a Displacement Discontinuity Method (DDM) based model. The results show that the sliding of fracture surface is irreversible but may change significantly after fluid pressure dissipates. The final sliding distance is different for natural and hydraulic fractures. Most of the shear fractures are natural fractures while the newly formed hydraulic fractures tend to be totally closed after pressure dissipates. The effects of in situ stress are investigated. The affected area reaches its maximum value when the maximum principle stress direction is perpendicular to the principal fracture direction. The effects of the injection rate are also investigated. The increasing of the injection rate is helpful in increasing the fracture aperture, but has no effect on the final sliding distance. Moreover, the effects of the injection rate on the affected area depend on the connectivity of natural fractures. The affected area increases with the injection rate when the connectivity is poor but decreases slightly with injection rate when the connectivity is good.


Introduction
Several mechanisms can lead to the permeability enhancement of the shale gas reservoir during hydraulic fracturing treatment [1]: (1) the propagation of Hydraulic Fractures (HFs); (2) the opening of Natural Fractures (NFs).This is reversible, i.e., the fractures will close after fluid pressure dissipates; and (3) the shear stimulation, and is irreversible.Clearly, the formation of a shear fracture network is important to the hydraulic fracturing treatment of shale gas reservoir.
However, the investigation on the formation of a shear fracture network is challenging.The hydraulic fracturing process depends on many important factors, such as natural fractures, rock properties, in situ stress, fluid properties, injection rate, etc.Moreover, the physical phenomena is not only complex but also ill conditioned [2].The simulation domain is often hundreds of meters, whereas the typical fracture aperture is a small fraction of a millimeter.A deformation that is considered small "noise" in the solid solver may induce dramatic oscillation of fluid pressure in the flow solver [2].Many models in early stages, such as Perkins-Kern (PK) model [3], Perkins-Kern-Nordgren (PKN) model [4], Khristianovich-Geertsma-Deklerk (KGD) mode [5], pseudo-3D models, and planar-3D models [6], cannot be used in simulating the propagation of complex fracture network.To meet the need of hydraulic fracturing design, many numerical models, such as the finite element method [7], extended finite element method [8], discrete element method [2] and mesh less method [9][10][11][12][13][14] have been developed in recent years.For example, Zhuang and Rabczuk et al. [12][13][14] proposed a meshless method that is able to treat the nucleation of fractures and complex patterns involving fracture Energies 2016, 9, 299 2 of 16 branching and crossing in a simple and effective way.Fu, Johnson and Carrigan [2] developed a model based on a discrete element method and simulated the propagation of fracture network with tens of fractures.The Displacement Discontinuity Method (DDM) [15] is another novel method for modelling the ill conditioned fluid-rock coupling system of hydraulic process.First, fracture displacements can be calculated with very high precision because the analytical solution is directly used to calculate the induced stress.Second, the grid number is much less than that in other methods because the rock matrix is not discretized.The propagation of complex fracture network with hundreds of fractures can be easily simulated by DDM based models [16][17][18][19][20][21][22].Moreover, based on the traditional DDM, Verde and Ghassemi [23,24] proposed a fast multipole displacement discontinuity method by considering the difference when calculating the induced stress of far field elements and the near field elements.Using this technique, a fracture system with up to 100,000 boundary elements can be simulated on currently mainstream computers.In summary, there are a variety of methods that can be used in modelling the hydraulic fracturing process nowadays.
Many researchers have investigated the shearing of fractures.For example, Pine and Batchelor [25] found that microseismic events are associated with shearing and the shearing during fluid injection is caused by the anisotropic in situ stress interacting with critically aligned joints.Willis-Richards et al. [26] presented a model in which the amount of shear displacement depends on the fracture shear stiffness and on the amount of "excess" shear stress available.Zhang et al. [27] derived the governing equations and the boundary conditions for equilibrium shear fractures, simulated the plane-strain fluid-driven fracture propagates and showed that the shear strength of a fracture can be reduced as fluid pressure increases until shearing is possible.De Bremaecker and Ferris [28] simulated the fracture propagation and showed that shear fractures can propagate in three different patterns, depending on the lateral normal stress.The fracture propagation occurs by a purely closed shear fracture at sufficiently high lateral normal stress.Majer et al. [29] showed that microseismic monitoring can only detect shear failure because of the very high frequency signals of tensile failure in many instances.Chipperfield et al. [30] presented a shear dilation diagnostics tool for evaluating tight gas stimulation treatments using the results of a simulation model which was devised to honor shear failure mechanisms.Nagel and Sanchez-Nagel [31] evaluated the shear failure along fracture surfaces as a function of fracture-induced stress and stress shadowing by utilizing the discrete element model.Nagel et al. [32] built a three-dimensional model of a reservoir section using a real Discrete Fracture Network (DFN), simulated the propagation of a large number of fractures and showed that the lower injection rate and the lower viscosity favor the creation of shear fracture.Zoback et al. [33] implemented both laboratory experiments and numerical modelling, and demonstrated the slow slip on pre-existing fractures is important to the effectiveness of slick-water hydraulic fracturing in shale gas reservoirs.Jung [34] reviewed the results and observations of the major Enhanced Geothermal System (EGS) projects and proved that the basic mechanism controlling the EGS projects is not the shearing of the fracture network but the formation of single large wing-cracks.Zangeneh et al. [35] investigated the relationship between hydraulic fracturing and the triggering of fracture slip of a naturally fractured rock mass.Riahi and Damjanac [1] performed a series of comparative studies and found that higher injection rates result in smaller shear stimulated areas.McClure and Horne [36] investigated the conditions required for shear stimulation and showed that hydraulic fracturing stimulation may occur through a mixture of the opening and sliding of preexisting fractures and the propagation of new fractures.McClure and Horne [20] proposed a tendency-for-shear-stimulation test for determining the degree to which shear stimulation contributes to stimulation and argued that shear stimulation is the only possible mechanism when bottom hole fluid pressure is slightly less than the minimum principal stress.In summary, the existing works have demonstrated the importance of the fracture shearing and have reached many useful conclusions.For example, the microseismic events are associated with shearing source mechanism [25], lower injection rate and the lower viscosity favor the creation of shear fracture [32], etc.Most of these conclusions are based on the fracture configurations during hydraulic fracturing treatment.However, the fracture geometries during fluid injection and shale gas production processes may not be the same because the fractures that have been created or reopened during fluid injection may partially close during production after fluid pressure dissipates [37].Therefore, it would be important to investigate the configurations of shear fracture network after fluid pressure dissipation.
In this paper, a DDM based model is developed and complex fracture networks are used to investigate the formation of shear fracture network.Considering that the fracture configurations may change after fluid pressure dissipates [37], both the hydraulic fracturing and the fluid leak off processes are simulated to get the final configurations of shear fracture network.We will first check how the fracture sliding distance changes after pressure dissipation and then investigate the optimization of the final shear fracture network after fluid pressure dissipation.

Numerical Method
A DDM-based model is used to simulate the hydraulic fracturing process.The following assumptions are used: the domain of rock matrix is infinite, and the rock matrix is homogeneous, isotropous and linear elastic [15].The rock matrix is treated as impermeable due to the ultra-low permeability of shale matrix.The fluid injected is Newtonian, single phase and laminar [20,38,39].
The elastic equilibrium equations for a system of N fracture elements are [39]: 0 rG 11 px, sq w psq `G12 px, sq ν psqs K px, sq ds where x = (x, y) is the coordinate, w is the normal displacement discontinuity, v is the shear displacement discontinuity, l r is the length of fracture r, and G ij is the hyper singular Green's functions [39].σ n is the normal stress and τ s is the shear stress, obeying Coulomb's frictional law characterized by the coefficient of friction λ: which can act in parts of fractures that are in contact, and vanishes along the separated parts.Along the opened fracture portions with fluid pressure, p f , we have: K is the three dimensional correction coefficient.Using the parameters given by Wu and Olson [40], the correction coefficient K proposed by Olson [41] can be written as: where h is the limite layer thickness perpendicular to the simulation plain, and d is the distance between points x 1 and x 2 .Volumetric flux q of fluid in fracture is given by: q " ´pw `w0 q 3 12µ where µ is fluid dynamic viscosity, w 0 , which is a reflection of the fracture surface roughness, is the initial hydraulic aperture, w is the fracture aperture, and Bs is the fluid pressure gradient.
Energies 2016, 9, 299 4 of 16 The global mass balance requires: where v f is fluid volume, L is fracture segment length, and 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 [39,42]: where K I and K II are stress intensity factors, K IC is tensile mode fracture toughness, θ is the fracture propagation direction relative to the current fracture orientation and satisfies [39,42]: 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.

Model Validation
The solving of the fracture aperture is one of the most important foundations for the modelling.Given the fluid pressure, the aperture of a single fracture in infinite rock matrix follows [15]: where a is the fracture half length, p f is the net fluid pressure, x is the distance to the center of the fracture, G is the shear modulus, and ν Po is 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 1.The numerical results agree well with the analytical one, even when there are only five grids.
The precise solving of the tip element aperture is especially important because the intensity factors are estimated by the displacements of the tip element [43]: where w and v are the normal and shear displacements of the tip element respectively, κ " 3 ´4v Po , and 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 use of the constant displacement element.However, the error is almost a constant.Therefore, a correction coefficient is introduced, i.e.,: where C = 0.79 is the correction coefficient.We refer for comparison to the model of Wu and Olson [44] for the fracture configurations of a single horizontal wellbore with two initial fractures (Figure 2a) and two horizontal wellbores with one fractures each (Figure 2b).These are the most classic fracture configurations and have been simulated in many works [18,19,32,40,44,45].The correct simulation of these configurations requires the precise 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 as that of Wu and Olson [44].The results are shown in Figure 3.For the initial configuration, as shown in Figure 3a, the two fractures propagate away from each other, while for the initial configuration as shown in Figure 3b, the two fractures grow towards each other.In summary, our model compares favorably with previous works.We refer for comparison to the model of Wu and Olson [44] for the fracture configurations of a single horizontal wellbore with two initial fractures (Figure 2a) and two horizontal wellbores with one fractures each (Figure 2b).These are the most classic fracture configurations and have been simulated in many works [18,19,32,40,44,45].The correct simulation of these configurations requires the precise 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 as that of Wu and Olson [44].The results are shown in Figure 3.For the initial configuration, as shown in Figure 3a, the two fractures propagate away from each other, while for the initial configuration as shown in Figure 3b, the two fractures grow towards each other.In summary, our model compares favorably with previous works.We refer for comparison to the model of Wu and Olson [44] for the fracture configurations of a single horizontal wellbore with two initial fractures (Figure 2a) and two horizontal wellbores with one fractures each (Figure 2b).These are the most classic fracture configurations and have been simulated in many works [18,19,32,40,44,45].The correct simulation of these configurations requires the precise 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 as that of Wu and Olson [44].The results are shown in Figure 3.For the initial configuration, as shown in Figure 3a, the two fractures propagate away from each other, while for the initial configuration as shown in Figure 3b, the two fractures grow towards each other.In summary, our model compares favorably with previous works.

Dimensionless Numbers
The anisotropy of crustal stress is important to the propagation of hydraulic fractures.Considering that the rock strength can be quantified by the initiation fluid pressure   , the stress anisotropy will be normalized by   in this paper.Here,   refers to the net fluid pressure above which HFs begin to propagate.As discussed by Zhang et al. [17],   is proportional to   √ ⁄ , where s is the average fracture spacing;, then, we define: where SA is the dimensionless stress anisotropy, σ max and σ min are the two principle stresses, respectively, and s is the average fracturing spacing of natural fracture network.
Specifically, when the principle stress directions are parallel to the coordinate axis, the dimensionless stress difference SA can be written as: where σ  and σ  are the principle stresses along x and y directions, respectively.The viscosity 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. [17]: where q is the injection rate, μ is fluid dynamic viscosity, and E is the Young's modulus.

Dimensionless Numbers
The anisotropy of crustal stress is important to the propagation of hydraulic fractures.Considering that the rock strength can be quantified by the initiation fluid pressure p ini , the stress anisotropy will be normalized by p ini in this paper.Here, p ini refers to the net fluid pressure above which HFs begin to propagate.As discussed by Zhang et al. [17], p ini is proportional to K IC { ?s, where s is the average fracture spacing;, then, we define: where SA is the dimensionless stress anisotropy, σ max and σ min are the two principle stresses, respectively, and s is the average fracturing spacing of natural fracture network.Specifically, when the principle stress directions are parallel to the coordinate axis, the dimensionless stress difference SA can be written as: where σ x and σ y are the principle stresses along x and y directions, respectively.The viscosity 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. [17]: where q is the injection rate, µ is fluid dynamic viscosity, and E is the Young's modulus.

Fracture Sliding During Pressure Dissipates
The fluid injection into a single natural fracture as shown in the sub graph of Figure 4a is first simulated.Fluid is injected from the center of the natural fracture.K IC is set to be infinite so the fracture propagation is forbidden.All of the other parameters are listed in Table 4.The variation of the average fracture aperture and sliding distance are shown in Figure 4a.Here, the fracture aperture is the normal displacement w in Equation ( 1), and the sliding distance is the shear displacement v in Equation (1).During the leak off process, the fracture aperture decreases.By contrast, the sliding distance is constant.

Fracture Sliding During Pressure Dissipates
The fluid injection into a single natural fracture as shown in the sub graph of Figure 4a is first simulated.Fluid is injected from the center of the natural fracture.KIC is set to be infinite so the fracture propagation is forbidden.All of the other parameters are listed in Table 4.The variation of the average fracture aperture and sliding distance are shown in Figure 4a.Here, the fracture aperture is the normal displacement w in Equation ( 1), and the sliding distance is the shear displacement v in Equation (1).During the leak off process, the fracture aperture decreases.By contrast, the sliding distance is constant.
The fluid injection into two natural fractures as shown in the sub graph of Figure 4b is also simulated.The results are different with the single fracture case.The average sliding distance of the two fractures is not a constant but decreases after fluid pressure dissipates.However, the sliding distance is greater than zero when the fractures are closed.These results indicate that the sliding of fracture is irreversible during fracturing treatment.Moreover, the sliding distance may change after fluid pressure dissipates.The coefficient of friction 0.9

Numerical Setting
The numerical setting is shown in Figure 5.The natural fracture network is reconstructed by a discrete fracture network model [1].The intersection angle between the natural fractures and x-axis is 60°.Fluid is injection with constant rate from the center.The fluid injection stops when any hydraulic fracture reaches the region boundary, i.e., the red dashed box in Figure 5.The pressure dissipation process is simulated based on the results of hydraulic fracturing simulation.KIC is set to be 1.0 × 10 6 Pa•m 0.5 , and the other default parameters are listed in Table 4, which are used in the following sections if there is no special statement.The fluid injection into two natural fractures as shown in the sub graph of Figure 4b is also simulated.The results are different with the single fracture case.The average sliding distance of the two fractures is not a constant but decreases after fluid pressure dissipates.However, the sliding distance is greater than zero when the fractures are closed.These results indicate that the sliding of fracture is irreversible during fracturing treatment.Moreover, the sliding distance may change after fluid pressure dissipates.

Numerical Setting
The numerical setting is shown in Figure 5.The natural fracture network is reconstructed by a discrete fracture network model [1].The intersection angle between the natural fractures and x-axis is 60 ˝.Fluid is injection with constant rate from the center.The fluid injection stops when any hydraulic fracture reaches the region boundary, i.e., the red dashed box in Figure 5.The pressure dissipation process is simulated based on the results of hydraulic fracturing simulation.K IC is set to be 1.0 ˆ10 6 Pa¨m 0.5 , and the other default parameters are listed in Table 4, which are used in the following sections if there is no special statement.

Effects of Stress Anisotropy
The hydraulic fracturing process under different values of SA are simulated.The principle crustal stress angle ψ is set to zero, i.e., the principle stresses are parallel to the coordinate axis.The fracture networks after fluid injection are shown in Figure 6.When SAxy = 0, the hydraulic fractures propagate farther along the y-axis direction.This is caused by the fact that the intersection angle between natural fractures and the y-axis is smaller than that of the x-axis.When SAxy < 0, the principal natural fracture direction is parallel to the principle stress direction, a narrower fracturing region is formed and the fracture network is less complex.By contrast, when SAxy > 0, a very complex fracture network is formed.
The pressure dissipation process is then simulated and the final shear fracture networks are shown in Figure 7.The fracturing region after pressure dissipation is comparable with that after fluid injection.However, the final sliding distance of a fracture has little relation with the fracture aperture and fluid pressure during fluid injection, i.e., the simulation of fluid pressure dissipation is meaningful to the prediction of fracture network configuration in production.
In Figure 7, most of the shear fractures are natural fractures.By contrast, although there are many newly formed hydraulic fractures, both the sliding and the aperture of these newly formed fractures tend to be disappear after pressure dissipates.This is because hydraulic fractures will change propagation direction when the fracture surface slides, which is required by the maximum hoop stress criterion that is defined by Equations ( 7) and (8).
The variation of the affected area, which is defined as the area of the region that experienced fluid pressure increase due to the fluid injection [1], is shown in Figure 8a.The affected area increases with SAxy, which is consistent with the results shown in Figure 6.Moreover, with the increase of stress difference, as shown in Figure 8b,c, both the average fracture density and the sliding distance in the fracturing region increase, which is beneficial to shale gas production.
The hydraulic fracturing processes under different principle stress angle ψ are simulated.The stress anisotropy SA is set to 0.5.The results are shown in Figure 9.It is clear that the affected area decreases with ψ.Moreover, there is a minimum value of fracture density and sliding distance when ψ ≈ 60°, which is the angle of one set of natural fractures.This indicates that it is unfavorable for shear simulation when the maximum principle stress direction is parallel to the natural fractures, under which circumstance the shear stress induced by crustal stress along natural fractures equals zero.

Effects of Stress Anisotropy
The hydraulic fracturing process under different values of SA are simulated.The principle crustal stress angle ψ is set to zero, i.e., the principle stresses are parallel to the coordinate axis.The fracture networks after fluid injection are shown in Figure 6.When SA xy = 0, the hydraulic fractures propagate farther along the y-axis direction.This is caused by the fact that the intersection angle between natural fractures and the y-axis is smaller than that of the x-axis.When SA xy < 0, the principal natural fracture direction is parallel to the principle stress direction, a narrower fracturing region is formed and the fracture network is less complex.By contrast, when SA xy > 0, a very complex fracture network is formed.
The pressure dissipation process is then simulated and the final shear fracture networks are shown in Figure 7.The fracturing region after pressure dissipation is comparable with that after fluid injection.However, the final sliding distance of a fracture has little relation with the fracture aperture and fluid pressure during fluid injection, i.e., the simulation of fluid pressure dissipation is meaningful to the prediction of fracture network configuration in production.
In Figure 7, most of the shear fractures are natural fractures.By contrast, although there are many newly formed hydraulic fractures, both the sliding and the aperture of these newly formed fractures tend to be disappear after pressure dissipates.This is because hydraulic fractures will change propagation direction when the fracture surface slides, which is required by the maximum hoop stress criterion that is defined by Equations ( 7) and (8).
The variation of the affected area, which is defined as the area of the region that experienced fluid pressure increase due to the fluid injection [1], is shown in Figure 8a.The affected area increases with SA xy , which is consistent with the results shown in Figure 6.Moreover, with the increase of stress difference, as shown in Figure 8b,c, both the average fracture density and the sliding distance in the fracturing region increase, which is beneficial to shale gas production.
The hydraulic fracturing processes under different principle stress angle ψ are simulated.The stress anisotropy SA is set to 0.5.The results are shown in Figure 9.It is clear that the affected area decreases with ψ.Moreover, there is a minimum value of fracture density and sliding distance when ψ « 60 ˝, which is the angle of one set of natural fractures.This indicates that it is unfavorable for shear simulation when the maximum principle stress direction is parallel to the natural fractures, under which circumstance the shear stress induced by crustal stress along natural fractures equals zero.First, we simulated the fracturing processes under different values of injection rate with poorly connected natural fracture network.The principle crustal stress angle ψ is set to 3π{4.The final fracture network configurations are shown in Figure 10.It is clear that the principal hydraulic fractures are formed when the value of M is small.By contrast, the large value of M will induce the very complex fracture network.Therefore, the affected area increases with the increasing of M, as shown in Figure 11a.The main reason is the localization instability of hydraulic fracture propagation can be partly suppressed by fluid viscosity when the injection rate is high.The stability problem of hydraulic fracture propagation has been discussed in our previous work [17].
The fracture propagation is dominated by the directions of both the natural fracture and the maximum principle stress.It can also be seen in Figure 10.The hydraulic fractures propagate along the maximum principle stress direction.However, when a hydraulic fracture encounters a natural fracture, it may be arrested by the natural fracture and then propagate along the natural fracture.The maximum principle stress direction may be different with the natural fracture direction.Therefore, the propagation direction of fracture network is controlled by both of the two directions.
The variation of fracture density and fracture displacement discontinuities are shown in Figure 11b,c.It is clear that the increasing of the injection rate is helpful for the increasing of fracture density and the displacement discontinuities during fracturing treatment.However, the increasing of M has no effect on the increasing of final sliding distance after pressure dissipates.These results indicate that the sliding of fractures is mainly controlled by the stress condition and the direction of natural fractures.
The fracturing processes with well-connected natural fracture networks are also simulated.The final fracture network configurations are shown in Figure 12.With the increase of M, the fracture density increases.This could also be seen from Figure 13b.However, the affected area does not increase with M as the poorly connected natural fracture network case.By contrast, the affected area may decrease, as shown in Figure 13a.This can be explained by the final fracture network configuration.When M is very big, i.e., the viscosity effect is strong, the shape of the fracturing region tends to be a round region that is centered with the injection point.By contrast, when M is small, i.e., the viscosity effect is weak, the fracturing fluid has sufficient time to fill most of the connected natural fractures.Therefore, the affected area may be larger when the M is smaller.The variation of fracture density and average sliding distance of fracture surface is similar with that of the poorly connected natural fracture network case.Both the fracture aperture and the sliding distance in fracturing treatment increases with the increasing of M.However, the final sliding distance of fracture surface does not change with the value of M.    16.       16.

Conclusions
The sliding of fracture surface is irreversible during the hydraulic fracturing process.However, the shear fracture network changes significantly after fluid pressure dissipates, i.e., the fracture network in fracturing treatment is different with the fracture network in shale gas production.Therefore, in fracturing treatment design, we should pay attention to the variation of fracture network during the pressure dissipation process.
The anisotropy of crustal stress is crucially important to the hydraulic fracture propagation and the final fracture configuration.The affected area reaches its maximum value when the maximum principle stress direction is perpendicular to the principal fracture direction.
The effect of natural fracture can be concluded into two aspects.First, the propagation direction of hydraulic fracture is controlled by the principal natural fracture direction.Second, as the sliding distance of hydraulic fracture is smaller than that of the natural fracture, the conductance of natural fracture network is important to the conductance of the final fracture network.Therefore, the natural fracture network is very important to the fracturing treatment and the production of shale gas.
The increasing of injection rate is helpful to the fracturing treatment by increasing the fracture aperture and the sliding distance during the fracturing process.However, the increasing of the injection rate has no effect on the final fracture sliding distance after fluid pressure dissipates.
The effects of the injection rate on the affected area depend on the connectivity of natural fracture.The affected area increases with the injection rate when the connectivity of the natural fracture is poor.However, when the connectivity of the natural fracture is good, the increasing of the injection rate has no effect on the increasing of the affected area.By contrast, the affected area may be smaller when the injection rate is high.Therefore, the injection rate should be determined according to the connectivity of the natural fracture during the design of hydraulic fracturing treatment.

Figure 2 .
Figure 2. 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 1 .
Figure 1.The comparison of numerical modelling with analytical solution.

Table 1 . 1 Figure 1 .
Figure 1.The comparison of numerical modelling with analytical solution.

Figure 2 .
Figure 2. 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 2 .
Figure 2. Two initiation fractures in (a) a horizontal wellbore and (b) two horizontal wellbores.Blue lines represent horizontal wellbores, and red lines represent initial fractures.

Table 3 . 16 Table 3 . 35 Figure 3 .
Figure 3.Comparison of propagation paths for initiation fractures in (a) a horizontal wellbore and (b) two horizontal wellbores.The line width represents the fracture aperture.

Figure 3 .
Figure 3.Comparison of propagation paths for initiation fractures in (a) a horizontal wellbore and (b) two horizontal wellbores.The line width represents the fracture aperture.

Figure 4 .
Figure 4.The variation of the average sliding distance and the average fracture aperture with time for (a) a single fracture and (b) two intersected fractures.The numerical settings are shown in the sub figures.The maximum principle crustal stress direction is represented by the angle ψ.The natural fractures are represented by the black segments.

Figure 4 .
Figure 4.The variation of the average sliding distance and the average fracture aperture with time for (a) a single fracture and (b) two intersected fractures.The numerical settings are shown in the sub figures.The maximum principle crustal stress direction is represented by the angle ψ.The natural fractures are represented by the black segments.

Figure 5 .
Figure 5.The numerical setting.The natural fractures are represented by the blue lines.The fluid is injected from the center.The boundary of the simulation region is shown as the red dashed box.The maximum principle stress direction is represented by the angle ψ.

Figure 5 .
Figure 5.The numerical setting.The natural fractures are represented by the blue lines.The fluid is injected from the center.The boundary of the simulation region is shown as the red dashed box.The maximum principle stress direction is represented by the angle ψ.

Figure 6 .
Figure 6.The fracture network configurations after fluid injection under different values of crustal stress anisotropy.The fluid pressure is represented by color.The fracture apereture is represented by curve width.(a)   = −0.5;(b)   = 0; (c)   = 0.5.

Figure 7 .
Figure 7.The final shear fracture networks after fluid pressure dissipates under different values of crustal stress anisotropy.The gray segments are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a)   = −0.5;(b)   = 0; (c)   = 0.5.

Figure 6 . 16 Figure 6 .
Figure 6.The fracture network configurations after fluid injection under different values of crustal stress anisotropy.The fluid pressure is represented by color.The fracture apereture is represented by curve width.(a) SA xy " ´0.5;(b) SA xy " 0 ; (c) SA xy " 0.5.

Figure 7 .
Figure 7.The final shear fracture networks after fluid pressure dissipates under different values of crustal stress anisotropy.The gray segments are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a)   = −0.5;(b)   = 0; (c)   = 0.5.

Figure 7 .Figure 8 .
Figure 7.The final shear fracture networks after fluid pressure dissipates under different values of crustal stress anisotropy.The gray segments are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a) SA xy " ´0.5;(b) SA xy " 0 ; (c) SA xy " 0.5.

Figure 9 .
Figure 9.The variation of (a) the affected area; (b) the fracture density; and (c) the average aperture and sliding distance with the direction of maximum principle stress.Here, "Final" refers to the values after fluid pressure dissipation.

Figure 8 . 16 Figure 8 .
Figure 8.The variation of (a) the affected area; (b) the fracture density; and (c) the average aperture and sliding distance with stress anisotropy SA xy .Here, "Final" refers to the values after fluid pressure dissipation.

Figure 9 .
Figure 9.The variation of (a) the affected area; (b) the fracture density; and (c) the average aperture and sliding distance with the direction of maximum principle stress.Here, "Final" refers to the values after fluid pressure dissipation.

Figure 9 .
Figure 9.The variation of (a) the affected area; (b) the fracture density; and (c) the average aperture and sliding distance with the direction of maximum principle stress.Here, "Final" refers to the values after fluid pressure dissipation.

Figure 10 .
Figure 10.The final fracture networks when the connectivity of natural fractures is poor.The gray segments are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a) M = 0.058; (b) M = 1.8;(c) M = 58.Here M is the dimensionless fluid viscosity that defined in Equation 16.

Figure 11 .
Figure 11.The variation of (a) the affected area; (b) the fracture density; and (c) the averaged displacements with the rate of injection when connectivity of natural fracture is poor.Here, "Final" refers to the values after fluid pressure dissipation.

Figure 10 .
Figure 10.The final fracture networks when the connectivity of natural fractures is poor.The gray segments are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a) M = 0.058; (b) M = 1.8;(c) M = 58.Here M is the dimensionless fluid viscosity that defined in Equation 16.

Energies 2016, 9 , 299 12 of 16 Figure 10 .
Figure 10.The final fracture networks when the connectivity of natural fractures is poor.The gray segments are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a) M = 0.058; (b) M = 1.8;(c) M = 58.Here M is the dimensionless fluid viscosity that defined in Equation 16.

Figure 11 .
Figure 11.The variation of (a) the affected area; (b) the fracture density; and (c) the averaged displacements with the rate of injection when connectivity of natural fracture is poor.Here, "Final" refers to the values after fluid pressure dissipation.Figure 11.The variation of (a) the affected area; (b) the fracture density; and (c) the averaged displacements with the rate of injection when connectivity of natural fracture is poor.Here, "Final" refers to the values after fluid pressure dissipation.

Figure 11 .
Figure 11.The variation of (a) the affected area; (b) the fracture density; and (c) the averaged displacements with the rate of injection when connectivity of natural fracture is poor.Here, "Final" refers to the values after fluid pressure dissipation.Figure 11.The variation of (a) the affected area; (b) the fracture density; and (c) the averaged displacements with the rate of injection when connectivity of natural fracture is poor.Here, "Final" refers to the values after fluid pressure dissipation.

Figure 12 .
Figure 12.The final fracture networks when the connectivity of natural fractures is good.The gray lines are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a) M = 0.058; (b) M = 1.8;(c) M = 58.Here M is the dimensionless fluid viscosity that defined in Equation 16.

Figure 13 .
Figure 13.The variation of (a) the affected area; (b) the fracture density; and (c) the averaged displacements with the rate of injection when connectivity of natural fracture is good.Here, "Final" refers to the values after fluid pressure dissipation.

Figure 12 .
Figure 12.The final fracture networks when the connectivity of natural fractures is good.The gray lines are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a) M = 0.058; (b) M = 1.8;(c) M = 58.Here M is the dimensionless fluid viscosity that defined in Equation 16.

Figure 12 .
Figure 12.The final fracture networks when the connectivity of natural fractures is good.The gray lines are natural fractures.The blue curves are hydraulic fractures.The shear fractures are colored red.(a) M = 0.058; (b) M = 1.8;(c) M = 58.Here M is the dimensionless fluid viscosity that defined in Equation 16.

Figure 13 .
Figure 13.The variation of (a) the affected area; (b) the fracture density; and (c) the averaged displacements with the rate of injection when connectivity of natural fracture is good.Here, "Final" refers to the values after fluid pressure dissipation.

Figure 13 .
Figure 13.The variation of (a) the affected area; (b) the fracture density; and (c) the averaged displacements with the rate of injection when connectivity of natural fracture is good.Here, "Final" refers to the values after fluid pressure dissipation.

Table 2 .
The comparison of numerical modelling with analytical solution for the aperture of tip element.Grid Number Numerical (m) Analytical (m) Analytical/Numerical

Table 2 .
The comparison of numerical modelling with analytical solution for the aperture of tip element.

Table 2 .
The comparison of numerical modelling with analytical solution for the aperture of tip element.