Comparison of Different Hydraulic Fracturing Scenarios in Horizontal Wells Using XFEM Based on the Cohesive Zone Method

Multistage hydraulic fracturing is a highly effective method for creating multiple transverse fractures to improve gas and oil reservoir production. It is critical to minimize the fracture spacing while also ensuring transverse propagation of fractures in multi-fractured horizontal wells. In this paper, a 3D fully coupled pore pressure-stress model based on the extended finite element method (XFEM) combined with the cohesive zone method is established to simulate five different fracturing scenarios in close spacing. The sensitivity of mesh size and the integration method are optimal, which are verified by the highly accurate traditional cohesive zone method. Then, the effect of five different fracturing scenarios on fracture geometries is compared. It is shown that spacing is a key parameter controlling fracture geometries in all fracturing scenarios. Alternative sequential and modified two-step fracturing can significantly reduce the influence of stress shadowing to generate more transverse fractures and form longer effective fractures. The sequential and two-step fracturing see an obvious improvement with increased fracture effective length when the spacing increases. The simultaneous fracturing technique can result in excessive closure of the middle fractures, which causes serious insertion of proppants. These results offer a new insight on optimization of hydraulic fracturing and can be a guidance for typical field cases.


Introduction
Unconventional reservoirs have ultralow permeability and porosity and often the stimulated reservoir volumes need to be increased to achieve commercial production rates.Multistage hydraulic fracturing in horizontal wells is a highly effective treatment for enhancing permeability, which places multiple transverse fractures along a horizontal well.When designing the fracturing scenario, it is important to minimize the spacing between clusters for higher production [1,2].However, closely-spaced clusters may have excessive interactions between fractures, which may result in undesirable fracture propagation patterns [3].The mechanical interaction between multiple fractures arises from a fracture's "stress shadowing" effect, which is mainly determined by the net pressure and the opening of propped fractures.Hence, great attention should be paid on the optimized trajectories of transverse fractures and to avoid excessive interference of stress.However, different hydraulic fracturing scenarios would redistribute the stress field and generate a variety of fracture networks.The mechanical mechanism needs a quantitative description with different hydraulic fracturing scenarios, in order to avoid undesirable diversion of fractures.
Due to the complex propagation patterns of fracture networks, it is difficult to measure the geometry of fractures in field applications using diagnostic techniques.Therefore, numerical simulations have become an important tool to solve this complex multiple field coupling problem.The problems associated with modeling multistage hydraulic fractures in horizontal wells have been addressed by many papers in recent years.By using an analytical solution, Solima et al. studied the effect of multiple fractured-induced interference, which would result in an increase in fracturing pressure [4].Wu and Olson developed a displacement discontinuity method (DDM) to simulate the propagation of multiple, potential, non-transverse fractures [5].This model couples the fracture deformation with fluid flow and can present the interaction effect between multiple fractures.Zou et al. simulated the arbitrary propagation of hydraulic fractures in naturally fractured formations, using the discrete element method (DEM) [6].In order to identify the plastic zone and the softening effect on fracture tips, Haddad et al. developed a finite element method (FEM) based on the cohesive zone method (CZM) to investigate the stress interference between horizontal clusters [7].The extend finite element method (XFEM) had been proposed to simulate the fracture propagation along arbitrary paths, which can be presented without the limitation of element re-meshing [8,9].This method adds the special shape function to the standard finite approximation with the framework of partition of unity.By adding "edge" phantom nodes with pore pressure degrees of freedom, the so alled XFEM based on CZM can be extended to simulate the hydraulic fracturing problems [10].Compared to the recent models for simulating multistage hydraulic fracturing in horizontal wells, this method provides the following improvements [11]: (1) fully coupled pore pressure-stress analysis and the effect of fractured-induced changes can be studied; (2) the leak-off behaviors are incorporated based on the fluid continuity and coupled with matrix deformations; and (3) more plausible modeling of quasi-brittle and heterogeneous reservoirs.
In the process of multistage hydraulic fracturing, the in-situ stress and cluster spacing play a vital role in creating fracture networks [12].Therefore, most studies focus on the optimization of the cluster spacing, which can avoid the stress interference from multiple fractures.Haddad et al. demonstrated the superiority of simultaneous fracturing to achieve the best fracture geometry by using the XFEM model.The study verified the strong interference of fractured-induced stress in a closely-spacing model [10].Liu et al. developed a XFEM model to investigate the sequential and simultaneous fracturing in a horizontal well, and the optimal spacing to create a large fracture network field was provided [13].Wang et al. studied the different propagation patterns of sequential and alternative sequential fracturing in horizontal wells [1].The XFEM model presents the stress orientation caused by the non-transverse fractures and describes the opening behaviors of propped hydraulic fractures.All the studies described above show that the complex interactions between fractures and different fracturing scenarios can result in different geometry of fractures.However, it is also a challenge for practical field applications without a clear understanding and optimization of different fracturing scenarios.
In this paper, a three dimensional (3D) model is developed to simulate nonplanar propagation of fractures in the porous media, which uses the XFEM method based on the CZM.Both pore-elasticity and seepage of rock matrix are taken into consideration.This model also allows us to simulate the leak-off between fracture surface and deformations.The mechanical interaction of propped fractures is simulated with five different types of fracturing scenarios: (1) simultaneous fracturing, (2) sequential fracturing, (3) alternative sequential fracturing, (4) two-step fracturing, and (5) modified two-step fracturing.The fracture network is characterized in a low in-situ stress scenario and uses close spacing.The optimal scenario is achieved by calculating the whole effective length of fracture network and comparing the geometries of each fracture network.

Fluid-solid Coupling Equation of Rock
In this simulation, the characteristic of seepage-stress is fully coupled, in which the fluid diffusion in porous medium is taken into account.Only a single phase of fluid is saturated in the pores of rock.The total stress can be determined by the effective stress by the rock skeleton and the pore pressure.According to the principle of virtual work, the equilibrium equation can be written as follows [14]: where the σ and p w are the effective stress and pore pressure, respectively; δε and δv are the virtual strain state and virtual velocity of deformation, respectively; t and f are the surface traction per unit area and body force per unit volume, respectively; I is the unit matrix.Fluid flow is simulated in the porous medium according to the continuity equation of liquid flow.The seepage behavior satisfies the typical Darcy's law.Therefore, the total change of fluid mass equals to the fluid mass crossing the porous medium per unit time.The continuity equation can be described as follows [15]: where the J is the change ratio of porous media; ρ w , n w and v w are the mass density of liquid, porosity of medium and average velocity of liquid relative to the solid phase, respectively.x is the space vector.Darcy's law presents the relation between the velocity of seepage flow and the gradient pressure in porous medium, which is written in the following form [16]: where k and g represent the hydraulic conductivity tensor and the gravity acceleration vector, respectively.

XFEM Based on Cohesive Zone Method
In order to avoid re-meshing during the fracture evolution, the XFEM method is introduced to model the arbitrary discontinuities in meshes, which results in a great improvement in fracture tip solutions.Compared to the traditional finite element method, the local enrichment functions are incorporated into a finite approximation function [8,9].Therefore, the extension of new displacement field u can be discretized into three parts based on the equilibrium equation.As shown in Figure 1, three parts of functions describe the nodes behaviors of whole elements, elements which are totally crossed and elements with crack-tip, respectively.Further, a fully coupled pore fluid diffusion and stress analysis has been implemented by adding pore pressure node.A complete presentation of this XFEM equation can be given in: where Sall is the set of all ordinary nodes, Sfrac and Stip are the set of Heaviside enrichment nodes and the set of fracture-tip enrichment nodes, respectively.u I is node degrees of freedom in regular finite element.a I and b l I represent the enrichment nodal degrees of freedom for fracture and crack-tip, 1, where x is the signed distance function which determines the two sides of a crack.
The third term is used to detail the nodes which are only cut by the crack tip.The asymptotic crack-tip functions Fl(x) is applied to describe the displacement around the fracture tip [17]: where r and θ are the polar coordinates at the crack tip with its origin coordinate system.

Crack Initiation and Damage Evolution
Unlike the methods used in previous studies [18][19][20][21][22], the CZM method is applied because of its potential superiority in modeling softening and plastic effects of fracture tips.The constitutive behavior of the CZM is defined by the traction-separation law, which includes three steps of facture propagation: (1) loading initiation; (2) damage initiation; and (3) damage evolution.
Damage initiation occurs at the beginning of stiffness degradation of the cohesive material.The process of degradation begins when the strains and stress satisfy the critical value in the assigned criteria.In this simulation, the maximum principle is chosen: where 0 max is the maximum allowable principal stress.The symbol〈〉represents the Macaulay bracket, which is used to indicate a pure compressive stress does not initiate damage.It assumes that damage will initiate when the stress ratio f reaches the value of 1.
Once the damage initiates, a scalar variable, D, is introduced to represent the average damage in the process of fracture propagation.D evolves from 0 to 1 after damage initiation and the element stiffness has a fully degradation when D reaches a value of one.Three components of one normal stress and two shear stress are determined linearly by the damage as follows [23]: (1 ) , 0 , 0 (1 ) where n t , s t and t t are the stress component in traction separation behavior without damage.
Here, the index n, s, t represent the direction of normal stress and the two shear stress, respectively.The second term on the right-hand of Equation ( 4) represents the nodes which are cut by the crack interior.These nodes are called Heaviside enriched nodes and the function H(x) is called the jump shape function: where x is the signed distance function which determines the two sides of a crack.The third term is used to detail the nodes which are only cut by the crack tip.The asymptotic crack-tip functions F l (x) is applied to describe the displacement around the fracture tip [17]: where r and θ are the polar coordinates at the crack tip with its origin coordinate system.

Crack Initiation and Damage Evolution
Unlike the methods used in previous studies [18][19][20][21][22], the CZM method is applied because of its potential superiority in modeling softening and plastic effects of fracture tips.The constitutive behavior of the CZM is defined by the traction-separation law, which includes three steps of facture propagation: (1) loading initiation; (2) damage initiation; and (3) damage evolution.
Damage initiation occurs at the beginning of stiffness degradation of the cohesive material.The process of degradation begins when the strains and stress satisfy the critical value in the assigned criteria.In this simulation, the maximum principle is chosen: where σ 0 max is the maximum allowable principal stress.The symbol represents the Macaulay bracket, which is used to indicate a pure compressive stress does not initiate damage.It assumes that damage will initiate when the stress ratio f reaches the value of 1.
Once the damage initiates, a scalar variable, D, is introduced to represent the average damage in the process of fracture propagation.D evolves from 0 to 1 after damage initiation and the element stiffness has a fully degradation when D reaches a value of one.Three components of one normal stress and two shear stress are determined linearly by the damage as follows [23]: t t = (1 − D)t t (10) where t n , t s and t t are the stress component in traction separation behavior without damage.Here, the index n, s, t represent the direction of normal stress and the two shear stress, respectively.Many formations are considered as quasi-brittle materials, in which linear softening effect is more suitable to simulate the process zone of fracture tip.Therefore, the variable, D, can be written linearly as follows: where δ m is the relative effective displacement attained over the course of loading.The superscripts max, f and 0 represent the max relative displacement, displacement when fractures are completely damaged, and the displacement at damage initiation.
To account for the displacement of three direction of stress, the effective separation is defined as: where δ n, δ s and δ t are the displacement in the direction of one normal stress and two shear stresses, respectively.Damage evolution can be described according to the fracture energy theory, when the damage occurs.The Benzeggagh-Kenane (BK) criterion is applied in this simulation to model the fracture propagation.Figure 2 shows fracture evolution used by the mix-mode BK criterion.The corresponding formulas is expressed as [24]: where G represents the energy release rate.Superscript c is the critical energy release rate.The subscripts n, s, t denote the directions of the normal stress and the two shear stresses, and η is a constant for a corresponding material property.Many formations are considered as quasi-brittle materials, in which linear softening effect is more suitable to simulate the process zone of fracture tip.Therefore, the variable, D, can be written linearly as follows: where δm is the relative effective displacement attained over the course of loading.The superscripts max, f and 0 represent the max relative displacement, displacement when fractures are completely damaged, and the displacement at damage initiation.
To account for the displacement of three direction of stress, the effective separation is defined as: where δn, δs and δt are the displacement in the direction of one normal stress and two shear stresses, respectively.Damage evolution can be described according to the fracture energy theory, when the damage occurs.The Benzeggagh-Kenane (BK) criterion is applied in this simulation to model the fracture propagation.Figure 2 shows fracture evolution used by the mix-mode BK criterion.The corresponding formulas is expressed as [24]: where G represents the energy release rate.Superscript c is the critical energy release rate.The subscripts n, s, t denote the directions of the normal stress and the two shear stresses, and η is a constant for a corresponding material property.

Fluid Flow within CZM
As shown in Figure 3, both normal and tangential flow are taken into consideration in the CZM fractures.The normal flow presents the crossing behavior between fracture surface and rock deformation.The tangential flow is the fluid flow along the fracture, which can describe the gradient pressure along the fracture path.All fluid within the model are assumed to be Newtown fluid which is incompressible and governed by the lubrication equation [25].

Fluid Flow within CZM
As shown in Figure 3, both normal and tangential flow are taken into consideration in the CZM fractures.The normal flow presents the crossing behavior between fracture surface and rock deformation.The tangential flow is the fluid flow along the fracture, which can describe the gradient pressure along the fracture path.All fluid within the model are assumed to be Newtown fluid which is incompressible and governed by the lubrication equation [25].The volume balance equation of tangential flow can be described by the following expression: where q is the tangential flow rate, w is the fracture width, μ is the viscosity of the fracturing fluid, and ∇pf is the gradient of pressure along the cohesive zone.This equation assumes that the injection rate is proportional to the pressure gradient and the cube of the fracture width and inversely proportional to the fracturing fluid viscosity.The normal flow is described as: where qt and qb are the fluid rate of permeation for the up and side surfaces, respectively; ct and cb are the leak-off coefficients of the top and bottom cohesive layers, respectively; pt and pb are the pore pressures on the top and bottom surfaces, respectively; pi is the pore pressure in the middle of the cohesive element.Therefore, the conservation equation of fluid mass containing normal and tangential flow is summarized as [26]: where Q is the total volume of injection fluid.

Model Construction and Verification
Hydraulic fracturing using XFEM based on the CZM method is more expensive in computation, when deviation of fractures occurs.This model needs a very fine mesh in a large field of enrichment which satisfies a requirement of quick convergence of solutions [27].All the simulations are solved by the standard solver based on the ABAQUS.In order to optimize the computing resource, we model a 3D hydraulic fracturing model using a single layer of fine element.Further, the sensitivity of mesh size and accuracy of integration method are investigated comparing with the FEM based on the CZM model, which has high accuracy in modeling hydraulic fracturing [25,[28][29][30].All parameters used in simulations come from the oil field of East China, as shown in Table 1.
ABAQUS presents a 3D solid element with pore pressure to model 3D hydraulic fracturing problems, which are solved by the analysis steps of geostatic and soils.Two main integration methods, including full integration (C3D8P) and reduced integration (C3D8RP) are provided to solve the problems coupling multiple factors.However, the sensitivity of mesh size and the accuracy of integration method were not verified in previous studies [3,11,13,31].In this paper, we use a single layer whose thickness is 0.5 m.Then, the size of horizontal is refined for four types as shown in Figure 4. Further, each size of C3D8P and C3D8RP are simulated to compare with the FEM based on CZM, which has a high accuracy in modeling hydraulic fractures [15,30,[32][33][34].Figures 5 and 6 show that the mesh size influences the accuracy of simulation results obviously.As the refining of elements, the injection pressure and mouth opening displacement (MOD) are closer to the results of FEM based on The volume balance equation of tangential flow can be described by the following expression: where q is the tangential flow rate, w is the fracture width, µ is the viscosity of the fracturing fluid, and ∇p f is the gradient of pressure along the cohesive zone.This equation assumes that the injection rate is proportional to the pressure gradient and the cube of the fracture width and inversely proportional to the fracturing fluid viscosity.The normal flow is described as: where q t and q b are the fluid rate of permeation for the up and side surfaces, respectively; c t and c b are the leak-off coefficients of the top and bottom cohesive layers, respectively; p t and p b are the pore pressures on the top and bottom surfaces, respectively; p i is the pore pressure in the middle of the cohesive element.Therefore, the conservation equation of fluid mass containing normal and tangential flow is summarized as [26]: where Q is the total volume of injection fluid.

Model Construction and Verification
Hydraulic fracturing using XFEM based on the CZM method is more expensive in computation, when deviation of fractures occurs.This model needs a very fine mesh in a large field of enrichment which satisfies a requirement of quick convergence of solutions [27].All the simulations are solved by the standard solver based on the ABAQUS.In order to optimize the computing resource, we model a 3D hydraulic fracturing model using a single layer of fine element.Further, the sensitivity of mesh size and accuracy of integration method are investigated comparing with the FEM based on the CZM model, which has high accuracy in modeling hydraulic fracturing [25,[28][29][30].All parameters used in simulations come from the oil field of East China, as shown in Table 1.
ABAQUS presents a 3D solid element with pore pressure to model 3D hydraulic fracturing problems, which are solved by the analysis steps of geostatic and soils.Two main integration methods, including full integration (C3D8P) and reduced integration (C3D8RP) are provided to solve the problems coupling multiple factors.However, the sensitivity of mesh size and the accuracy of integration method were not verified in previous studies [3,11,13,31].In this paper, we use a single layer whose thickness is 0.5 m.Then, the size of horizontal is refined for four types as shown in Figure 4. Further, each size of C3D8P and C3D8RP are simulated to compare with the FEM based on CZM, which has a high accuracy in modeling hydraulic fractures [15,30,[32][33][34].Figures 5 and 6 show that the mesh size influences the accuracy of simulation results obviously.As the refining of elements, the injection pressure and mouth opening displacement (MOD) are closer to the results of FEM based on CZM.The horizontal size refines to 0.5 × 0.5 m can achieve high accuracy both using the C3D8P and the C3D8RP.Owing to high computation speed of reduced integration, several scholars applied the C3D8RP in hydraulic fracturing simulations [35,36].CZM.The horizontal size refines to 0.5 × 0.5 m can achieve high accuracy both using the C3D8P and the C3D8RP.Owing to high computation speed of reduced integration, several scholars applied the C3D8RP in hydraulic fracturing simulations [35,36].Nevertheless, the injection pressure and mouth opening displacement using the C3D8P are more accurate comparing with the results from the C3D8RP, as shown in Figures 5 and 6.It indicates that the results of C3D8P have a higher accuracy than the results of C3D8RP with the requirement of mesh refinement.Further, we discovered that C3D8RP has a disadvantage in solution convergence which will result in tremendous distortion of elements, as shown in Figure 7.
Refining the mesh size can reduce the convergence failures but increasing the mesh numbers, which will increase the computational cost.Therefore, we restrict the element for C3D8P and the Nevertheless, the injection pressure and mouth opening displacement using the C3D8P are more accurate comparing with the results from the C3D8RP, as shown in Figures 5 and 6.It indicates that the results of C3D8P have a higher accuracy than the results of C3D8RP with the requirement of mesh refinement.Further, we discovered that C3D8RP has a disadvantage in solution convergence which will result in tremendous distortion of elements, as shown in Figure 7.
Refining the mesh size can reduce the convergence failures but increasing the mesh numbers, which will increase the computational cost.Therefore, we restrict the element for C3D8P and the optimal mesh size is 0.5 × 0.5 × 0.5 m in this work.In order to keep the injection of fluid when the fracture initiates, the initial gap of enrichment is set at 2 mm.In order to reduce the computing cost, only one half of the model is simulated by using a symmetric model.As shown in Figure 8, the computational domain has a length of 300 m and a width of 200 m.All elements are set to be enrichment, but only 100 × 100 m domain adopts the optimal mesh size of 0.5 × 0.5 m to guarantee the accuracy of fracture propagation.In order to eliminate the effect of boundaries, the model extends far enough.Simultaneously, the domains without fracture  In order to reduce the computing cost, only one half of the model is simulated by using a symmetric model.As shown in Figure 8, the computational domain has a length of 300 m and a width of 200 m.All elements are set to be enrichment, but only 100 × 100 m domain adopts the optimal mesh size of 0.5 × 0.5 m to guarantee the accuracy of fracture propagation.In order to eliminate the effect of boundaries, the model extends far enough.Simultaneously, the domains without fracture  In order to reduce the computing cost, only one half of the model is simulated by using a symmetric model.As shown in Figure 8, the computational domain has a length of 300 m and a width of 200 m.All elements are set to be enrichment, but only 100 × 100 m domain adopts the optimal mesh size of 0.5 × 0.5 m to guarantee the accuracy of fracture propagation.In order to eliminate the effect of boundaries, the model extends far enough.Simultaneously, the domains without fracture In order to reduce the computing cost, only one half of the model is simulated by using a symmetric model.As shown in Figure 8, the computational domain has a length of 300 m and a width of 200 m.All elements are set to be enrichment, but only 100 × 100 m domain adopts the optimal mesh size of 0.5 × 0.5 m to guarantee the accuracy of fracture propagation.In order to eliminate the effect of boundaries, the model extends far enough.Simultaneously, the domains without fracture propagation used the transition mesh sizes, which can reduce the computational Energies 2019, 12, 1232 9 of 19 cost effectively.Therefore, the entire domain is discretized into 54830 elements.A zero-displacement degree of freedom and an initial pore pressure (10 MPa) boundary are applied on the outer boundary (AB, BC, CD).Further, a symmetric boundary is set at the horizontal well (AD) and the displacement of all nodes is fixed in the thickness-direction.In order to compare the effect of different fracturing scenarios, three fracturing stages are adopted in this simulation.The numbers 1 , 2 , and 3 represent the right, middle and left fracture, respectively.
Energies 2018, 11, x FOR PEER REVIEW 9 of 19 propagation used the transition mesh sizes, which can reduce the computational cost effectively.Therefore, the entire domain is discretized into 54830 elements.A zero-displacement degree of freedom and an initial pore pressure (10 MPa) boundary are applied on the outer boundary (AB, BC, CD).Further, a symmetric boundary is set at the horizontal well (AD) and the displacement of all nodes is fixed in the thickness-direction.In order to compare the effect of different fracturing scenarios, three fracturing stages are adopted in this simulation.The numbers ① , ② , and ③ represent the right, middle and left fracture, respectively.

Simultaneous Fracturing
The case of simultaneous fracturing in horizontal wells is investigated in this section, in which the fracture spacings are 10 and 20 m.The used parameters are from an oilfield in East China, as shown in Table 1.The geometries for simultaneous fracturing fractures are shown in Figure 9.The propagation of multiple fractures can influence the fracture geometries in a variety of patterns from coalescence to divergence.It indicates that the stress field changes a lot in the near region of fracture surface and fracture-tip.This phenomenon, which is induced by the propped fracture and the diffusion of pore-pressure, can result in undesirable propagation of fracture networks and reduces the efficiency of each fracturing stage.
Unlike the common assumption in previous studies that the trajectories of fractures will propagate in the direction of maximum horizontal principle stress [25,37], with simultaneous fracturing, only the middle fracture follows a direction that is perpendicular to the minimum horizontal principle stress.The side fractures tend to deviate, which is due to the stress interference induced by the middle fracture.The side fractures are suppressed to be short and wide compared with the middle fracture.Furthermore, the middle fracture can break though rapidly and negatively impact the propagation in the whole process of fracturing.Conversely, the fracture width of the middle fracture is influenced by the stress shadowing, which leads to excessive closure of the middle fracture in the near well region.The phenomenon produces an excessive pressure on the fracture surface and results in a narrow flow channel for the injection fluid.The final trajectories of fractures indicate that the side fractures deviate from the desired orthogonal paths and middle fracture suffered an excessive pressure near the wellbore.In this situation, the effectiveness of the three fractures decreases significantly.
Different results of injection pressure and MOD for the spacing at 10 and 20 m can be achieved (as shown in Figure 10a,b).From Figure 10a, we know that the injection pressure of side fractures are much higher than the middle fracture, which is caused by the stress shadowing.However, the injection pressure of the two side fractures are consistent.This indicates that the effect of side fractures can reduce the propagation pressure of the middle fracture, which is the main reason that

Simultaneous Fracturing
The case of simultaneous fracturing in horizontal wells is investigated in this section, in which the fracture spacings are 10 20 m.The used parameters are from an oilfield in East China, as shown in Table 1.The geometries for simultaneous fracturing fractures are shown in Figure 9.The propagation of multiple fractures can influence the fracture geometries in a variety of patterns from coalescence to divergence.It indicates that the stress field changes a lot in the near region of fracture surface and fracture-tip.This phenomenon, which is induced by the propped fracture and the diffusion of pore-pressure, can result in undesirable propagation of fracture networks and reduces the efficiency of each fracturing stage.With the spacing = 20 m, even the mouth of the middle fracture will close, and the average facture width is much higher than the case of spacing = 10 m.That can be explained by the fact that an increment of spacing can reduce the stress shadowing significantly in a closely-spaced scenario.When the spacing is increased, a longstanding opening of the fracture mouth can be maintained, which is useful to the injection and placement of proppant for simultaneous fracturing.Unlike the common assumption in previous studies that the trajectories of fractures will propagate in the direction of maximum horizontal principle stress [25,37], with simultaneous fracturing, only the middle fracture follows a direction that is perpendicular to the minimum horizontal principle stress.The side fractures tend to deviate, which is due to the stress interference induced by the middle fracture.The side fractures are suppressed to be short and wide compared with the middle fracture.Furthermore, the middle fracture can break though rapidly and negatively impact the propagation in the whole process of fracturing.Conversely, the fracture width of the middle fracture is influenced by the stress shadowing, which leads to excessive closure of the middle fracture in the near well region.The phenomenon produces an excessive pressure on the fracture surface and results in a narrow flow channel for the injection fluid.The final trajectories of fractures indicate that the side fractures deviate from the desired orthogonal paths and middle fracture suffered an excessive pressure near the wellbore.In this situation, the effectiveness of the three fractures decreases significantly.
Different results of injection pressure and MOD for the spacing at 10 and 20 m can be achieved (as shown in Figure 10a,b).From Figure 10a, we know that the injection pressure of side fractures are much higher than the middle fracture, which is caused by the stress shadowing.However, the injection pressure of the two side fractures are consistent.This indicates that the effect of side fractures can reduce the propagation pressure of the middle fracture, which is the main reason that generates a deeper propagation of the middle fracture.The fracture geometries of opening and closure can influence the efficiency of hydraulic fractures.Excessive closure of fractures will result in great difficult in fluid injection and increase the possibility of proppant insertion.In this situation, the drainage of the reservoirs decreases significantly.Figure 10b presents the variation of MOD in the fracturing process.Compared with side fractures, the MOD of middle fractures decreases significantly as the two side fractures evolve.The difference is that the maximum MOD in larger spacing (20 m) is wider than with spacing = 10 m.Further, in the spacing = 10 m case, the middle fracture opens when the fractures initiate from the well, and suffers excessive suppression in the early stages of fracturing.In that situation, the mouth of the middle fracture will completely close in a short period of time.With the spacing = 20 m, even the mouth of the middle fracture will close, and the average facture width is much higher than the case of spacing = 10 m.That can be explained by the fact that an increment of spacing can reduce the stress shadowing significantly in a closely-spaced scenario.When the spacing is increased, a longstanding opening of the fracture mouth can be maintained, which is useful to the injection and placement of proppant for simultaneous fracturing.generates a deeper propagation of the middle fracture.The fracture geometries of opening and closure can influence the efficiency of hydraulic fractures.Excessive closure of fractures will result in great difficult in fluid injection and increase the possibility of proppant insertion.In this situation, the drainage of the reservoirs decreases significantly.Figure 10b presents the variation of MOD in the fracturing process.Compared with side fractures, the MOD of middle fractures decreases significantly as the two side fractures evolve.The difference is that the maximum MOD in larger spacing (20 m) is wider than with spacing = 10 m.Further, in the spacing = 10 m case, the middle fracture opens when the fractures initiate from the well, and suffers excessive suppression in the early stages of fracturing.In that situation, the mouth of the middle fracture will completely close in a short period of time.With the spacing = 20 m, even the mouth of the middle fracture will close, and the average facture width is much higher than the case of spacing = 10 m.That can be explained by the fact that an increment of spacing can reduce the stress shadowing significantly in a closely-spaced scenario.When the spacing is increased, a longstanding opening of the fracture mouth can be maintained, which is useful to the injection and placement of proppant for simultaneous fracturing.

Sequential and Alternative Sequential Fracturing
In this section, the sequential and alternative sequential fracturing are simulated using different close spacings of 10 and 20 m.All of parameters are listed in Table 1.Two groups of fracturing scenarios are taken into account to optimize a desirable sequential fracturing treatment in a close spacing scenario.The final fracture geometries of sequential and alternative sequential fracturing are presented in Figure 11.

Sequential and Alternative Sequential Fracturing
In this section, the sequential and alternative sequential fracturing are simulated using different close spacings of 10 and 20 m.All of parameters are listed in Table 1.Two groups of fracturing scenarios are taken into account to optimize a desirable sequential fracturing treatment in a close spacing scenario.The final fracture geometries of sequential and alternative sequential fracturing are presented in Figure 11.
the propagation of the middle fracture (third fracture).Further, less deviation of fractures can be observed in the alternative sequential fracturing comparing with the sequential fracturing.This denotes that fracturing sequence can easily change the fracture geometries in a close spacing situation.In the case 2 (Figure 11b) and case 4 (Figure 11d), the geometries change a lot both in sequential and alternative sequential fracturing, when the spacing increases from 10 to 20 m.It denotes that spacing is a key parameter in controlling the fracture geometries.However, the left fracture propagates longer and the middle fracture presents more curvature, when the alternative sequential fracturing is applied.This demonstrates that the alternative sequence can enhance the interference on the middle fracture similar to the simultaneous fracturing, but alleviates the suppression on the left fracture.In order to investigate the efficiency of different fracturing scenarios, the variation of injection pressure and MOD with sequential and alternative sequential fracturing are shown in Figure 12.As shown in Figure 12a,c, all stages are fracturing continuously without a flow back of injection fluid, which can simulate the propping of proppants.Therefore, the injection pressure of each fracture In case 1 (Figure 11a), the sequential fracturing is applied in the situation of spacing = 10 m, whose fracturing sequence is 1 -2 -3 (right fracture-middle fracture-left fracture).In this situation, only the right fracture (first fracture) propagates along an orthogonal path and extends deepest compared with the subsequent fractured fractures.The middle fracture (second fracture) and left fractures (third fracture) deviate from the desirable direction significantly, and present obvious coalescence with the right fracture.Moreover, the middle fracture propagates deeper than the left fracture, which can be explained by the fact that the left fracture has suffered stronger influence from the prior fractured right and middle fracture.However, the middle fracture is only influenced by the right fracture.Owing to stress shadowing in the third fracturing stage, the right and middle fractures come to be closed near the wellbore.Simultaneously, the width of the right fracture decreases the most and tends to close completely in the near well region.The fracture geometries of case 3 (Figure 11c) are much different in case 1.In case 3, the alternative fracturing uses a scenario of spacing = 10 m, whose fracturing sequence is 1 -3 -2 (right fracture-left fracture-middle fracture).The alternative sequence can enhance the propagation of the left fracture (second fracture), but suppresses the propagation of the middle fracture (third fracture).Further, less deviation of fractures can be observed in the alternative sequential fracturing comparing with the sequential fracturing.This denotes that fracturing sequence can easily change the fracture geometries in a close spacing situation.In the case 2 (Figure 11b) and case 4 (Figure 11d), the geometries change a lot both in sequential and alternative sequential fracturing, when the spacing increases from 10 to 20 m.It denotes that spacing is a key parameter in controlling the fracture geometries.However, the left fracture propagates longer and the middle fracture presents more curvature, when the alternative sequential fracturing is applied.This demonstrates that the alternative sequence can enhance the interference on the middle fracture similar to the simultaneous fracturing, but alleviates the suppression on the left fracture.In order to investigate the efficiency of different fracturing scenarios, the variation of injection pressure and MOD with sequential and alternative sequential fracturing are shown in Figure 12.As shown in Figure 12a,c, all stages are fracturing continuously without a flow back of injection fluid, which can simulate the propping of proppants.Therefore, the injection pressure of each fracture remains constant, when the fracturing stage is completed (injection rate = 0 m 3 /s).When the initial fracture damages completely, the injection pressure of all cases decreases to a stable value which can guarantee the propagation of fractures.It is obvious that the latter fracturing stage has a higher injection pressure than previous fracturing stages, whether use the sequential or alternative sequential fracturing scenarios.Compared with the sequential fracturing, the second fractures in alternative fracturing have a lower injection pressure (close to 19.3 MPa), which is almost identical to the injection pressure of the first fractures (close to 18.8 MPa).However, the injection pressure of second fractures (close to 19.6 MPa) in sequential fracturing are closer to that of the third fractures (close to 20.2 MPa).Notably, the injection pressure of each fracturing stage is consistent when the spacing increases from 10 to 20 m.It demonstrates that the spacing hardly influences the injection pressure in a close spacing scenario, when using the same fracturing scenario.
pressure in a close spacing scenario, when using the same fracturing scenario.
Figure 12b,d show the MOD with sequential and alternative sequential fracturing.Without the influence of stress shadowing, the MOD curves show no difference between two fracturing scenarios in the first fracturing stage.The MOD of sequential fracturing fracture decreases faster than the alternative sequential fracturing in the latter fracturing stage.In this process, the MOD in the case of spacing = 20 m is always larger than the case of spacing = 10 m, which is due to the stronger influence of fracture-induced stress with shorter spacing.In the second and third fracturing stage, the MOD of first fractures decreases due to strong stress shadowing.In the situation of spacing = 10 m, the fracture mouth closes completely in the second fracturing stage (fracturing time of closure = 3850 s), which should be avoided for serious insertion of proppants, while in the situation of spacing =20 m, the first fracture remains open even all three fracturing stages are completed.The second fractured fractures experience a propagation pattern that is similar to the first fracture.However, a higher value of MOD can be achieved in the second fracturing stage, which is mainly caused by a higher value of the injection pressure.In this situation, the second fracture of sequential fracturing also has a larger MOD comparing with the alternative sequential fracturing.In the third fracturing stage of sequential fracturing, the second fracture in situation of spacing = 10 closes obviously.While in the situation of spacing = 20 m, the MOD of second fracture decreases in a short time and remains at a constant (close to 14.6 mm) value during the fracturing process.This may ensure a high efficient hydraulic fracture and reduces the possibility of proppant insertion.For the third fractures, all cases have a sharp increase of MOD except the case of sequential fracturing at spacing = 20 m.That can be explain that less stress interference leads to less consumption of injection rates in broadening the fracture.In that case, deeper propagation can be achieved for increasing the reservoir production.

Two-step Fracturing and Modified Two-step Fracturing
In order to enhance production in the reservoirs, it is important to minimize the fracture spacing while ensures the propagation of transverse fractures.In this section, the two-step and modified twostep fracturing are simulated to investigate the effect on the facture geometries, as shown in Figure Figure 12b,d show the MOD with sequential and alternative sequential fracturing.Without the influence of stress shadowing, the MOD curves show no difference between two fracturing scenarios in the first fracturing stage.The MOD of sequential fracturing fracture decreases faster than the alternative sequential fracturing in the latter fracturing stage.In this process, the MOD in the case of spacing = 20 m is always larger than the case of spacing = 10 m, which is due to the stronger influence of fracture-induced stress with shorter spacing.In the second and third fracturing stage, the MOD of first fractures decreases due to strong stress shadowing.In the situation of spacing = 10 m, the fracture mouth closes completely in the second fracturing stage (fracturing time of closure = 3850 s), which should be avoided for serious insertion of proppants, while in the situation of spacing = 20 m, the first fracture remains open even all three fracturing stages are completed.The second fractured fractures experience a propagation pattern that is similar to the first fracture.However, a higher value of MOD can be achieved in the second fracturing stage, which is mainly caused by a higher value of the injection pressure.In this situation, the second fracture of sequential fracturing also has a larger MOD comparing with the alternative sequential fracturing.In the third fracturing stage of sequential fracturing, the second fracture in situation of spacing = 10 closes obviously.While in the situation of spacing = 20 m, the MOD of second fracture decreases in a short time and remains at a constant (close to 14.6 mm) value during the fracturing process.This may ensure a high efficient hydraulic fracture and reduces the possibility of proppant insertion.For the third fractures, all cases have a sharp increase of MOD except the case of sequential fracturing at spacing = 20 m.That can be explain that less stress interference leads to less consumption of injection rates in broadening the fracture.In that case, deeper propagation can be achieved for increasing the reservoir production.

Two-step Fracturing and Modified Two-step Fracturing
In order to enhance production in the reservoirs, it is important to minimize the fracture spacing while ensures the propagation of transverse fractures.In this section, the two-step and modified two-step fracturing are simulated to investigate the effect on the facture geometries, as shown in Figure 13.In our simulation, two closely-spacing of 10, 20 m are used in the different two-step scenarios.The variation of injection pressure and MOD are calculated in this simulation as shown in Figure 14.In the first fracturing stage, the injection pressure decreases to a stable value (close to 19.4 MPa) after the initiation of first fractures and few differences are observed when two different fracturing scenarios are used (Figure 14a,c).In the second fracturing stage, the second fractures have a higher  In case 1 (Figure 13a), the spacing equals to 10 m and the fracturing sequence is 1 -2 2 (middle fracture-left and right fracture), which is called two-step fracturing.In this situation, the middle fracture (first fracture) initiates and propagates along a desirable orthogonal path.Then the side fractures are influenced by the middle fracture and tend to deviate seriously.As the fracturing progresses, the fracture tip of the side fractures reorient to the surface of the middle fracture, which will result in a more curved and shorter propagation path.Further, the fracture-induced stress of side fractures can conversely lead to excessive pressure on the surface of the middle fracture.Therefore, the middle fracture suffers obvious closure near the propagation path of the two side fractures.The similar propagation patterns can be observed in the condition of spacing = 20 m, but less deviation occurs in the larger spacing, as shown in Figure 13b.
In case 3 (Figure 13c), the modified two-step fracturing is applied in the condition of spacing = 10 m.The fracturing sequence is 1 1 -2 (right and left fracture-middle fracture).In this situation, the left and right fracture initiate and deviate from each other in the first fracturing stage.The lengths of the two side fractures are bigger than the length of the fracture in two-step fracturing and even an obvious stress interference effect can be observed.For the second fractured middle fracture, the growth is suppressed significantly in the whole process and the length is much smaller than the two side fractures.An undesirable effect of middle fracture is that the behavior of coalescence that leads an excessive closure for the right fracture, which will reduce the efficiency of the right fracture due to proppant insertion.With increasing spacing (spacing = 20 m (Figure 13d)), the less effect of stress interference realizes an average propagation of multiple fractures, especially for the middle fracture whose length is longer than in the condition of spacing = 10 m.
The variation of injection pressure and MOD are calculated in this simulation as shown in Figure 14.In the first fracturing stage, the injection pressure decreases to a stable value (close to 19.4 MPa) after the initiation of first fractures and few differences are observed when two different fracturing scenarios are used (Figure 14a,c).In the second fracturing stage, the second fractures have a higher injection pressure during the whole process compared with the first fracturing stage.For two-step fracturing, the injection pressure of second fractures is bigger than it is in modified two-step fracturing.Using different spacing = 10, 20 m, the injection pressure changes less.The two fracturing scenarios have different behaviors on the effect of MOD.When the middle fractures are first fractured (Figure 14b), the MOD of the middle fracture increases in the first fracturing stage, while it decreases to the value of the initial opening gap (2 mm) in the latter fracturing stage.During the process of decrease, the MOD in the condition of spacing = 10 m is smaller and reduces quicker than when spacing = 20 m.However, the MOD of the first fractures using modified two-step fracturing are quite different from the results of two-step fracturing, as shown in Figure 14d.The variation of injection pressure and MOD are calculated in this simulation as shown in Figure 14.In the first fracturing stage, the injection pressure decreases to a stable value (close to 19.4 MPa) after the initiation of first fractures and few differences are observed when two different fracturing scenarios are used (Figure 14a,c).In the second fracturing stage, the second fractures have a higher injection pressure during the whole process compared with the first fracturing stage.For two-step fracturing, the injection pressure of second fractures is bigger than it is in modified two-step fracturing.Using different spacing = 10, 20 m, the injection pressure changes less.The two fracturing scenarios have different behaviors on the effect of MOD.When the middle fractures are first fractured (Figure 14b), the MOD of the middle fracture increases in the first fracturing stage, while it decreases to the value of the initial opening gap (2 mm) in the latter fracturing stage.During the process of decrease, the MOD in the condition of spacing = 10 m is smaller and reduces quicker than when spacing = 20 m.However, the MOD of the first fractures using modified two-step fracturing are quite different from the results of two-step fracturing, as shown in Figure 14d.Both right and left fractures have different performance on MOD whether the spacing is 10 or 20 m.In the first fracturing stage of spacing = 10 m, we calculate that the MOD of left fracture experiences a rapidly increase and achieves a biggest value (16.18 mm), while the MOD of right fracture increases slowly and even starts to decrease in the later period of the first fracturing stage.That can be explained by the fact that the middle fracture has an obvious coalescence into the right fracture, which results in excessive suppression on the fracture wall of the right fracture.This phenomenon is observed in Figure 13c and should be avoided.In the second fracturing stage, all cases of MOD increase similarly and are consistent with each other.

Discussion
When considering the five fracturing scenarios together, the geometries of multiple fractures can be significantly influenced by different fracturing scenarios in a close spacing scenario.Different propagation patterns from fracture coalescence to divergence are formed because of fracture-induced stress shadowing.Using the simultaneous fracturing in close spacing, the side fractures experience obvious coalescence and outward deviation, while the middle fracture forms a transverse fracture.This phenomenon is similar to the results from Haddad et al. [10].Further, the middle fracture can break though and gets away from the influenced region of side fractures rapidly.Therefore, the middle fracture is much longer than the side fractures.The side fractures are suppressed seriously by the middle fracture and tend to form more curved fractures, which are shorter and wider.Inversely, the excessive opening of side fractures will result in excessive closure of the middle fracture along the propagation paths of the side fractures.
Compared with sequential fracturing, the alternative sequential fracturing can reduce the interference of stress shadowing and form more transverse fractures.Increasing the spacing will alleviate the influence both in sequential and alternative sequential fracturing, which is consistent with the conclusions stated by Wang et al. [1].Using an alternative fracturing, the optimal spacing can be reduced again.In this paper, we found that the later fracturing stage has a higher injection pressure compared with earlier fracturing stage, and the latter fractured fractures experience more serious suppression in the fracture growth.This indicates that increasing the number of fracturing stages is harmful to increasing the volume of the fracture network when close spacing is used.Further, the latter fractured fractures will result in excessive closure of previous fractures, which Both right and left fractures have different performance on MOD whether the spacing is 10 or 20 m.In the first fracturing stage of spacing = 10 m, we calculate that the MOD of left fracture experiences a rapidly increase and achieves a biggest value (16.18 mm), while the MOD of right fracture increases slowly and even starts to decrease in the later period of the first fracturing stage.That can be explained by the fact that the middle fracture has an obvious coalescence into the right fracture, which results in excessive suppression on the fracture wall of the right fracture.This phenomenon is observed in Figure 13c and should be avoided.In the second fracturing stage, all cases of MOD increase similarly and are consistent with each other.

Discussion
When considering the five fracturing scenarios together, the geometries of multiple fractures can be significantly influenced by different fracturing scenarios in a close spacing scenario.Different propagation patterns from fracture coalescence to divergence are formed because of fracture-induced stress shadowing.Using the simultaneous fracturing in close spacing, the side fractures experience obvious coalescence and outward deviation, while the middle fracture forms a transverse fracture.This phenomenon is similar to the results from Haddad et al. [10].Further, the middle fracture can break though and gets away from the influenced region of side fractures rapidly.Therefore, the middle fracture is much longer than the side fractures.The side fractures are suppressed seriously by the middle fracture and tend to form more curved fractures, which are shorter and wider.Inversely, the excessive opening of side fractures will result in excessive closure of the middle fracture along the propagation paths of the side fractures.
Compared with sequential fracturing, the alternative sequential fracturing can reduce the interference of stress shadowing and form more transverse fractures.Increasing the spacing will alleviate the influence both in sequential and alternative sequential fracturing, which is consistent with the conclusions stated by Wang et al. [1].Using an alternative fracturing, the optimal spacing can be reduced again.In this paper, we found that the later fracturing stage has a higher injection pressure compared with earlier fracturing stage, and the latter fractured fractures experience more serious suppression in the fracture growth.This indicates that increasing the number of fracturing stages is harmful to increasing the volume of the fracture network when close spacing is used.Further, the latter fractured fractures will result in excessive closure of previous fractures, which leads to Energies 2019, 12, 1232 16 of 19 a serious insertion of proppant.This undesirable effect can be alleviated by using the alternative fracturing treatment.
Kumar et al. [3] found that the two-step fracturing can be used in a relatively closely spaced scenario and enhance the reservoir permeability efficiently.The same conclusion was reached in our work.Using the two-step and modified two-step fracturing, the geometries of fractures can be obviously changed.In two-step fracturing, the middle fracture propagates along a desirable orthogonal path to form transverse fracture.While the side fractures have an obvious behavior of coalescence to the middle fracture, especially under the condition of close spacing = 10 m.However, the side fractures can propagate more transversely and longer when the modified two-step fracturing is used, but the middle fracture can be suppressed because of the two side fractures.In conclusion, the second fractured fractures always suffers an influence of the first fractured fractures.The modified fracturing produces an obvious improvement of the flow efficiency of fractures in a close spacing scenario.
In order to better understand the effect of different fracturing scenarios, the effective fracture length is provided in Figure 15.We consider a propagation depth of fractures as the effective fracture length.When the spacing increases, the effective fracture length increases in all fracturing scenarios.This indicates that spacing is a key factor in controlling the fracture geometries, even with close spacing.When the close spacing = 10 m, alternative sequential and modified two-step fracturing have an obvious advantage compared with other fracture scenarios, whose effective fracture length can reach 156 and 155 m, respectively.However, when the spacing increases to 20 m, the sequential and two-step fracturing show a large increase in effective fracture length, while the effective fracture length of modified two-step and alternative fracturing change less in this situation.Therefore, in a close spacing scenario (10 m), the alternative and modified fracturing scenario are recommended to generate more efficient fractures.However, the modified two-step can achieve less improvement of fracture length when the spacing increases to 20 m.This means that increasing spacing seems meaningless when the modified two-step fracturing is applied in a close spacing environment.
Energies 2018, 11, x FOR PEER REVIEW 16 of 19 leads to a serious insertion of proppant.This undesirable effect can be alleviated by using the alternative fracturing treatment.Kumar et al. [3] found that the two-step fracturing can be used in a relatively closely spaced scenario and enhance the reservoir permeability efficiently.The same conclusion was reached in our work.Using the two-step and modified two-step fracturing, the geometries of fractures can be obviously changed.In two-step fracturing, the middle fracture propagates along a desirable orthogonal path to form transverse fracture.While the side fractures have an obvious behavior of coalescence to the middle fracture, especially under the condition of close spacing = 10 m.However, the side fractures can propagate more transversely and longer when the modified two-step fracturing is used, but the middle fracture can be suppressed because of the two side fractures.In conclusion, the second fractured fractures always suffers an influence of the first fractured fractures.The modified fracturing produces an obvious improvement of the flow efficiency of fractures in a close spacing scenario.
In order to better understand the effect of different fracturing scenarios, the effective fracture length is provided in Figure 15.We consider a propagation depth of fractures as the effective fracture length.When the spacing increases, the effective fracture length increases in all fracturing scenarios.This indicates that spacing is a key factor in controlling the fracture geometries, even with close spacing.When the close spacing = 10 m, alternative sequential and modified two-step fracturing have an obvious advantage compared with other fracture scenarios, whose effective fracture length can reach 156 and 155 m, respectively.However, when the spacing increases to 20 m, the sequential and two-step fracturing show a large increase in effective fracture length, while the effective fracture length of modified two-step and alternative fracturing change less in this situation.Therefore, in a close spacing scenario (10 m), the alternative and modified fracturing scenario are recommended to generate more efficient fractures.However, the modified two-step can achieve less improvement of fracture length when the spacing increases to 20 m.This means that increasing spacing seems meaningless when the modified two-step fracturing is applied in a close spacing environment.

Conclusions
Multistage hydraulic fracturing in horizontal wells is developed based on a 3D-XFEM, which couples the pore pressure-stress combining the cohesive zone method.Five different fracturing scenarios, including simultaneous, sequential, alternative sequential, two-step and modified twostep fracturing are simulated to investigate the impact of fracture-induced interference on fracture geometries in close spacing scenarios.The results of the five different scenarios are compared.Based on the simulations, the following conclusions can be reached: (1) Fracture spacing is a key parameter in controlling the geometries of multiple fractures.Narrower spacing results in stronger interference of stress shadowing, which leads to a shorter

Conclusions
Multistage hydraulic fracturing in horizontal wells is developed based on a 3D-XFEM, which couples the pore pressure-stress combining the cohesive zone method.Five different fracturing scenarios, including simultaneous, sequential, alternative sequential, two-step and modified two-step fracturing are simulated to investigate the impact of fracture-induced interference on fracture geometries in close spacing scenarios.The results of the five different scenarios are compared.Based on the simulations, the following conclusions can be reached: (1) Fracture spacing is a key parameter in controlling the geometries of multiple fractures.Narrower spacing results in stronger interference of stress shadowing, which leads to a shorter fracture length and wider fracture width for latter fractured fractures.The latter fractured fractures also have an undesirable effect on excessive closure of previously fractured fractures.
(2) The stress shadowing in the propagation domain is redistributed when different scenarios are used.In a close spacing scenario (10 m), alternative sequential and modified two-step fracturing can significantly reduce the influence of stress shadowing to generate more transverse fractures and longer effective fracture lengths.
(3) The sequential and two-step fracturing have an obvious improvement in increasing fracture effective length when the spacing increases.Thus, more transverse and less closure of fractures will occur by increasing the spacing, when the two fracturing scenarios are used.
(4) The simultaneous fracturing technique can result in an excessive closure of the middle fracture, which causes serious insertion of proppants.In the later fracturing process, the injection of fluid will be difficult.The larger spacing should be optimal to reduce the influence on this situation.
The numerical model presented in this paper can be a guidance in optimizing the fracturing scenarios in closely-spacing horizontal wells.
u is the standard finite shape function of node I. H(x) is the enrichment shape jump function.F l is the singular displacement field around the facture tip.

Figure 1 .
Figure 1.Schematic of enrich nodes and common fluid nodes (pore-pressure nodes within a fracture).

Figure 1 .
Figure 1.Schematic of enrich nodes and common fluid nodes (pore-pressure nodes within a fracture).

Figure 5 .Figure 6 .
Figure 5. Variation of injection pressure and mouth opening displacement with C3D8P: (a) Variation of injection pressure; (b) Variation of mouth opening displacement.

Figure 5 .Figure 5 .Figure 6 .
Figure 5. Variation of injection pressure and mouth opening displacement with C3D8P: (a) Variation of injection pressure; (b) Variation of mouth opening displacement.

Figure 6 .Figure 5 .Figure 6 .
Figure 6.Variation of injection pressure and mouth opening displacement with C3D8RP: (a) Variation of injection pressure; (b) Variation of mouth opening displacement.

Figure 8 .
Figure 8. Construction of element model.

Figure 8 .
Figure 8. Construction of element model.
Energies 2018, 11, x FOR PEER REVIEW 10 of 19generates a deeper propagation of the middle fracture.The fracture geometries of opening and closure can influence the efficiency of hydraulic fractures.Excessive closure of fractures will result in great difficult in fluid injection and increase the possibility of proppant insertion.In this situation, the drainage of the reservoirs decreases significantly.Figure10bpresents the variation of MOD in the fracturing process.Compared with side fractures, the MOD of middle fractures decreases significantly as the two side fractures evolve.The difference is that the maximum MOD in larger spacing (20 m) is wider than with spacing = 10 m.Further, in the spacing = 10 m case, the middle fracture opens when the fractures initiate from the well, and suffers excessive suppression in the early stages of fracturing.In that situation, the mouth of the middle fracture will completely close in a short period of time.

Figure 10 .
Figure 10.Variation of injection pressure and mouth opening displacement with different closely-spacing simultaneous fracturing: (a) Variation of injection pressure; (b) Variation of mouth opening displacement.

Figure 11 .
Figure 11.Fracture geometries of sequential and alternative sequential fracturing with different closely-spacing (POR is the pore pressure: Pa): (a) Fracture geometries of sequential fracturing with spacing = 10 m; (b) Fracture geometries of sequential fracturing with spacing = 20 m; (c) Fracture geometries of alternative sequential fracturing with spacing = 10 m; (d) Fracture geometries of alternative sequential fracturing with spacing = 20 m.

Figure 11 .
Figure 11.Fracture geometries of sequential and alternative sequential fracturing with different closely-spacing (POR is the pore pressure: Pa): (a) Fracture geometries of sequential fracturing with spacing 10 m; (b) Fracture geometries of sequential fracturing with spacing = 20 m; (c) Fracture geometries of alternative sequential fracturing with spacing = 10 m; (d) Fracture geometries of alternative sequential fracturing with spacing = 20 m.

Figure 12 .
Figure 12.Variation of injection pressure and mouth opening displacement with different closely-spacing sequential and alternative sequential fracturing: (a) Variation of injection pressure with sequential fracturing; (b) Variation of mouth opening displacement with sequential fracturing; (c) Variation of injection pressure with alternative sequential fracturing; (d) Variation of mouth opening displacement with alternative sequential fracturing.

Energies 2018 , 19 Figure 13 .
Figure 13.Fracture geometries of two-step and modified two-step fracturing with different closelyspacing (POR is the pore pressure: Pa): (a) Fracture geometries of two-step fracturing with spacing = 10 m; (b) Fracture geometries of two-step fracturing with spacing = 20 m; (c) Fracture geometries of modified two-step fracturing with spacing = 10 m; (d) Fracture geometries of modified two-step fracturing with spacing = 20 m.

Figure 13 .
Figure 13.Fracture geometries of two-step and modified two-step fracturing with different closely-spacing (POR is the pore pressure: Pa): (a) Fracture geometries of two-step fracturing with spacing = 10 m; (b) Fracture geometries of two-step fracturing with spacing = 20 m; (c) Fracture geometries of modified two-step fracturing with spacing = 10 m; (d) Fracture geometries of modified two-step fracturing with spacing = 20 m.

Energies 2018 , 19 Figure 13 .
Figure 13.Fracture geometries of two-step and modified two-step fracturing with different closelyspacing (POR is the pore pressure: Pa): (a) Fracture geometries of two-step fracturing with spacing = 10 m; (b) Fracture geometries of two-step fracturing with spacing = 20 m; (c) Fracture geometries of modified two-step fracturing with spacing = 10 m; (d) Fracture geometries of modified two-step fracturing with spacing = 20 m.

Figure 14 .
Figure 14.Variation of injection pressure and mouth opening displacement with different closely-spacing two-step and modified two-step fracturing; (a) Variation of injection pressure with two-step fracturing; (b) Variation of mouth opening displacement with two-step fracturing; (c) Variation of injection pressure with modified two-step fracturing; (d) Variation of mouth opening displacement with modified two-step fracturing.

Figure 15 .
Figure 15.Variation of effective fracture length with different closely-spacing fracturing scenarios.

Figure 15 .
Figure 15.Variation of effective fracture length with different closely-spacing fracturing scenarios.

Table 1 .
Input parameters of simulations.

Table 1 .
Input parameters of simulations.