2016-04 Optimal spacing of sequential and simultaneous fracturing in horizontal well

Optimal spacing of sequential and simultaneous fracturing in horizontal wellChuang Liu a ,XiaoLong Wang a ,DaWei Deng b ,YongPing Zhang b ,YuGuang Zhang b ,HengAn Wu a ,*,He Liu c ,**aCAS Key Laboratory of Mechanical Behavior and Design of Materials,Department of Modern Mechanics,University of Science and Technology of China,Hefei,Anhui 230027,China bPetroChina Daqing Oil Field Limited Company,Daqing,Heilongjiang 163454,China cPetroChina Research Institute of Petroleum Exploration &Development,Beijing 100083,Chinaa r t i c l e i n f oArticle history:Received 10December 2015Received in revised form 4January 2016Accepted 18January 2016Available online 19January 2016Keywords:Unconventional reservoirs Optimal spacing Horizontal well Fracture network Sequential fracturing Simultaneous fracturinga b s t r a c tSequential and simultaneous fracturing are two main fracturing methods for creating multiple hydraulic fractures in horizontal wells.The fracture spacing of these two fracturing scenarios is critical for creating fracture network to improve unconventional gas or oil reservoirs production.In this paper,a fully coupled (flow and mechanics)hydraulic fracture propagation model based on extended finite element method (XFEM)is established to account for the sequential and simultaneous propagation of nonplanar fractures.By mapping the distribution of in-situ stress contrast during fracturing process,we calculate the fracture network area that is the extent of low in-situ stress contrast zones.The fracture spacing that satis fies the requirement of preventing sand plug as well as creating large field of fracture network is optimal.Sensitivity studies have been conducted for optimal spacing.It is shown that optimal spacing will decrease with in-situ stress contrast,and is not sensitive to varying Young's modulus of the rock matrix.The method presented in this paper can be used in horizontal well design to estimate reasonable fracture spacing.©2016Elsevier B.V.All rights reserved.1.IntroductionMultiple-fracture treatments in horizontal wells are the key technology that depletes unconventional oil or gas reservoirs economically (Pope et al.,2010;Chaudhri,2012;Yu and Sepehrnoori,2013).Due to the dif ficulties of diagnostic tech-niques for measuring individual fracture and fracture network,numerical simulation remains an important tool for engineers to estimate the geometry of fractures and fracture network.The problem associated with modeling hydraulic fractures has been addressed by many papers in recent years.Wu and Olson (2015)developed a numerical model based on displacement disconti-nuity method (DDM)to simulate simultaneous multi-fracture treatments for fully coupled fluid flow and fracture mechanics in horizontal wells.The cohesive element method has been applied to modeling viscosity-dominated hydraulic fractures (Zhang et al.,2010;Chen,2012;Wang et al.,2012,2015a ).Weng et al.(2014)presented a complex fracture network model that simulates hy-draulic fracture networks created during the stimulation treatment and proppant placement.Sesetty and Ghassemi (2015)developed displacement discontinuity method for sequential and simulta-neous hydraulic fracturing in single and multi-lateral horizontal wells.Wang (2015)proposed a numerical model for nonplanar hydraulic fracture propagation in brittle and ductile rocks using XFEM with cohesive zone method.The mathematical and numer-ical models for estimating the production of shale gas have been presented (Yu et al.,2014;Zhang et al.,2014,2015).There are many earlier works focus on creating fracture network.Olson and Taleghani (2009)studied multiple hydraulic fractures growth and their interaction with natural fractures based on displacement discontinuity method.They found that low in-situ stress anisotropy tends to improve fracture complexity in naturally fractured reser-voirs.Fu et al.(2013)built an explicitly coupled hydro-geomechanical model for simulating hydraulic fracturing in arbi-trary discrete fracture networks.Their simulation results demon-strate that large field of fracture network could be created in low in-*Corresponding author.**Corresponding author.E-mail addresses:wuha@ (H.Wu),liuhe@ (H.Liu).Contents lists available at ScienceDirectJournal of Natural Gas Science and Engineeringjournal h omepage:ww /locate/jngse/10.1016/j.jngse.2016.01.0241875-5100/©2016Elsevier B.V.All rights reserved.Journal of Natural Gas Science and Engineering 29(2016)329e 336situ stress contrast reservoirs.Guo et al.(2015)simulated a hy-draulic fracture crossing a cemented natural fracture.The in-situ stress contrast was found to be the key parameter that contains how a hydraulic fracture propagates across a natural fracture.The previous works showed the complex interaction between fractures and suggested that in-situ stress contrast plays a vital role in the creating of fracture network.However,the simulation of fracture network during fracturing process for practical engineering appli-cations is still a challenge work.The spacing between fractures is critical for the success of multiple fracture treatments of horizontal wells.There are many earlier works focusing on the interaction of fractures and spacing optimization.Roussel and Sharma(2011a)investigated the stress perturbation induced by the propped fracture.Their results demonstrated that a stress reorientation of90 can occur in the vicinity of a transverse fracture.This zone is called the stress reversal region,and the fracture spacing should be large enough to restrain the propagation of longitudinal fractures.They further investigated multiple-fracture deflection in horizontal well,and proposed that the spacing that makes the main fracture to deflect up to5 is optimal(Roussel and Sharma,2011b).Their results provided some new insights on spacing optimization in horizontal well.Morrill and Miskimins(2012)investigated the sensitivity of stress shadow to various reservoir properties and optimized hy-draulic fracture spacing in unconventional shales.Soliman et al. (2010)studied the Texas-two step method in horizontal well and presented that fracture spacing is designed to achieve isotropy of in-situ stress and facilitate the reorientation of in-situ stress.Yu et al.(2014)performed a sensitivity study of gas production for a shale gas well with different geometries of multiple transverse hydraulic fractures.They revealed that the outer fractures contribute more to gas production when fracture spacing is small due to the effect of fracture interference.Liu et al.(2015)numeri-cally simulated the connectivity of fracture network produced by Texas-two step method and proposed a new method for fracture spacing optimization in horizontal shale-gas well.The previous works focus on the interaction between fractures with a given fracture length.However,the in-situ stress is changing during the fracture propagation,which in turn affects the propagation and intersection of natural fractures and then the fracture network.The combination of nonplanar fractures propagation and fracture network during fracturing process for spacing optimization of sequential and simultaneous fracturing in horizontal well has not been presented yet.In the process of fracturing operation,the proppant is injected in conjunction with the fracturingfluid to prevent the closure of fracture surface after releasing thefluid pressure.The proppant is a number of small particles,typically sand,treated sand or man-made ceramic materials,so the fracture width must be large enough to make the proppant injection feasible.Zhou et al.(2014) investigated a low-efficient hydraulic fracturing operation in a tight gas reservoir in the North German Basin.The results demonstrated that the transporting and setting of proppant can result in a different closure of fracture surface.The full closure of a fracture at the perforation will lead to a low productivity.In practical,the fracture width has to be large enough along the propagation path when the proppant is added in the injectionfluid,otherwise sand plug will occur and results in serious engineering accident.It is necessary to investigate the mechanical interaction between frac-tures in horizontal well and then the width of the fractures to avoid the sand plug.In this paper,a numerical model of nonplanar fractures propa-gation in porous media was established to optimize fracture spacing in sequentially and simultaneously fractured horizontal well.The mechanical interaction of proppant to fracture surface and the sand plug during fracturing process are incorporated in our modeling and simulation.The fracture network is characterized by low in-situ stress contrast during fracture propagation.Optimal fracture spacing is achieved by calculating the average fracture network width with the varying fracture spacing in our model. 2.Modeling and simulation method2.1.Fracture propagation modelThe fracture propagation model presented in this paper is based on two-dimensional plane strain and couples porous media deformation andfluidflow.The extendedfinite element method (XFEM)is used to represent the mechanics of fractures and their opening,including the mechanical interaction with closely spaced fractures.Fluidflow in the fractures is determined by lubrication equation.The effect of propped fractures on in-situ stress as well as the propagation of other fractures is taken into account by applying truss model on fracture surfaces.The maximum in-situ compres-sive stress principle is used to determine the fracture propagation direction.2.1.1.The cohesive lowThe fracture initiation and propagation analysis is governed by cohesive traction-separation constitutive behavior.Linear elastic behavior followed by the initiation and evolution of damage is assumed in the traction-separation model.The elastic behavior is written in terms of an elastic constitutive matrix that relates the nominal stresses to the nominal strains across the interface.The nominal stresses are the force components divided by the original area,while the nominal strains are the separations divided by the original thickness.Damage initiation is the beginning of degradation of the response of a material.The process of degradation begins when the stresses or strains satisfy certain specified damage initiation criteria.The maximum principal stress criterion is adopted in this paper and it can be written asf¼h s max is0max(1)where,s0max represents the maximum allowable principal stress. The symbol〈〉represents the Macaulay bracket(i.e.,h s max i¼0if s max<0and h s max i¼s max if s max!0).Damage is assumed to initiate when the maximum principal stress ratio(as defined in the expression above)reaches a value of one.Damage evolution describes the degrading of material stiffness after the corresponding initiation criterion is reached.A scalar damage variable,D,represents the averaged overall damage.It initially has a value of0.If damage evolution occurs,D mono-tonically evolves from0to1after the initiation of damage.The stress components are affected by the damage according tot¼ð1ÀDÞt damage initiatedt no damage occurs(2)where t are stress components predicted by the linear elastic traction-separation behavior for the current separations without damage.For linear softening,the scalar damage variable,D,is expressed with the following formC.Liu et al./Journal of Natural Gas Science and Engineering29(2016)329e336 330D¼d f nd maxnÀd0nd maxnd f nÀd0n(3)where d f n and d0n are the opening displacements at complete failure and at the initiation of damage,respectively.d maxnrefers to the maximum value of opening displacement attained during the loading history.2.1.2.Fluidflow within a fractureMany fracturingfluid within fractures exhibit power law rheo-logical.In this paper,to enhance convergence stability with little cost in accuracy,incompressible and Newtonianfluid is assumed in this study.The tangentialflow rate q within the fracture is deter-mined by the lubrication equations(Batchelor,1967),which can be written asq¼Àw312mV p f(4)where p f is thefluid pressure inside the fracture,w is the fracture width,m is thefluid viscosity.The conservation offluid mass inside the fracture can be described asv wv tþV qþðq tþq bÞ¼0(5)where q t and q b are the normalflow rates into the two fracture surfaces,respectively.The normalflow is defined as8 <:q t¼c tp fÀp tq b¼c bp fÀp b(6)where p t and p b are pore pressures in the adjacent formation, respectively and c t and c b define the correspondingfluid leak-off coefficients.2.1.3.Fluidflow in porous mediaInfluidfilled porous media,the total stress s related to the ef-fect stress s0is expressed ass¼s0þa p w(7) where p w is pore pressure in porous media.a is poroelastic constant that is independent of thefluid properties.The a is assumed to be1 due to its negligible influence on fracture geometry(Ghassemi, 1996).Equilibrium equation can be written in the form of virtual work principal for the volume under its current configuration at time t asZ Vs0Àp w I:dεdV¼ZSt$d v dSþZVf$d v dV(8)where s0and dεare the effective stress and virtual rate of defor-mation respectively.t and f are the surface traction per unit area and body force per unit volume,I is unit matrix.Fluidflow in porous media satisfies the mass conservation equation written as (Malvern,1969)vðr w fÞv tþV$ðr w f v wÞ¼0(9)where r w is the density of porousfluid,f is the porosity and v w is the seepageflow velocity.According to Darcy's law,the relation between the velocity of seepageflow and the gradient of porous pressure is linear in the following form(Marino and Luthin,1982)v w¼À1f g r wk$ðp wÀr w gÞ(10)where k,p w and g represent the hydraulic conductivity tensor,the porous pressure and the gravity acceleration vector respectively.2.1.4.Extendedfinite element methodThe extendedfinite element method(XFEM)wasfirst intro-duced by Belytschko and Black(1999)to model arbitrary discon-tinuities in meshes.This extension exploits the partition of unity property offinite elements,which allows local enrichment func-tions to be easily incorporated into afinite element approximation while preserving the classical displacement variational setting (Melenk and Babuska,1996).The XFEM models the discontinuity ina displacementfield along the crack path,wherever this path maybe located with respect to the mesh.Thisflexibility enables the method to simulate crack growth without remeshing.The presenceof discontinuities is ensured by the special enriched functions in conjunction with additional degrees of freedom.We shall here summarize the main idea of this method.A complete presentationis given in detail(Dolbow and Belytschko,1999).In order to model the presence of the discontinuity,the displacement vector u can be approximated as(Moes and Belytschko,2002)u¼Xi2Iu i f iþXj2Jb j f j HðfðxÞÞþXk2Kf kXl¼1c lkF lðxÞ!(11)The second term on the right-hand side is valid for nodes whose shape function support is cut by the crack interior.The jump function Hð$Þis written asHðxÞ¼þ1x!0À1x<0(12)and fðxÞis the signed distance function to the crack(the sign determining whether x is on one side or the other of the crack).The third term is used only for nodes whose shape function support is cut by the crack tip.The asymptotic crack-tip functions F lðxÞare used to model the displacementfield around the tip of the discontinuity.These functions are chosen based on the asymptotic features of the displacementfield at the crack tipF lðxÞ¼ffiffiffirpsinq2;ffiffiffirpcosq2;ffiffiffirpsinq2sinðqÞ;ffiffiffirpcosq2sinðqÞ(13) whereðr;qÞis a polar coordinate system with its origin at the crack tip.Accurately modeling the crack-tip singularity requires constantly keeping track of where the crack propagates and is cumbersome because the degree of crack singularity depends on the location of the crack.Therefore,the asymptotic singularity functions is not considered in this model.2.2.Bulk modulus evolutionDuring fracturing process,the initiation and propagation of hydraulic fractures will change in-situ stress in its neighborhood. The variations of in-situ stress causes the changes of Young's modulus in the rock matrix according to poroelasticity theory. Zimmerman(1990)and Hassanzadegan and Zimmermann(2014) have expressed the drained bulk modulus for fractured porous media as an exponential function depending on the TerzaghiC.Liu et al./Journal of Natural Gas Science and Engineering29(2016)329e336331effective stress as1¼1∞þ1KÀ1∞eÀs e=p0(14)where K i and K∞refer to the bulk modulus at low and high effective pressures respectively.s e is the minimum effective in-situ stress.p0 is a characteristic closure pressure that depends on the solid and bulk properties.For isotropy matrix,Young's modulus can be determined by the bulk modulus asE¼3ð1À2vÞK(15) where K is calculated using Eq.(14).v is Poisson's ratio,and it is assumed to be constant due to its narrow range of variations.2.3.Truss modelIn hydraulic fracturing process,the proppant is injected in conjunction with the fracturingfluid to prevent the closure of the fracture after releasing thefluid pressure.The proppant is a great number of small particles which subjected to the confining stress and prevents the fracture from fully closing,leading to the residual opening of the fracture.The opening of a propped fracture in hor-izontal well causes a reorientation of in-situ stress in its neigh-borhood,which in turn affects the propagation of subsequent main fractures as what is called stress shadowing(Roussel and Sharma, 2011b;Taghichian et al.,2014).A linear elastic truss element model is used in this paper to simulate the mechanical interaction of settled proppant to fracture surface.The truss model is many truss elements distributed on opened fracture surface.The distributed truss element can experience closure on the application of external stresses.The closure of fracture surface is prevented by the truss elements due to the normal stress acting on it.The truss elements do not put on fracture surfaces initially.It is active after the fracture opened characterizing the injected proppant at the end of fracturing stage.2.4.Simulation of fracture networkSuccessfully creating fracture network in horizontal wells is critical for unconventional reservoirs production economically.In some Chinese oilfields,the formation in-situ stress contrast is large, which suppresses the growth of fracture network.In-situ stress contrast could be a controlling factor for creating fracture network in these reservoirs.During fracturing process,fluid loss to the for-mation will reopen cemented natural fractures in a zone where the stress contrast is under a critical value.The initiation and propa-gation of natural fractures are not considered in this paper.In-situ stress contrast is used to determine whether fracture network can be produced(if the in-situ stress contrast less than a given threshold value)in a region during stimulation process.The critical value of in-situ stress contrast is assumed to be3MPa in this paper.The average fracture network width(AFNW)is used in this paper.It is defined as fracture network area between fractures/ fracture spacing after fracturing.It can be used to characterize the regions of fracture network at different spacing.Fig.1illustrates the fracture network and the AFNW for two fractures in a horizontal well.The red(in the web version)region represents the produced fracture network between fractures after fracturing.3.Numerical results and discussion3.1.Model verificationA simulation case for model verification is presented in this section.The used parameters are from a Chinese oilfield.The in-jection rate is4.5m3/min with the period of300min.The hori-zontal in-situ stresses are50MPa and61MPa.The calculated surface pressure during fracturing operation were compared with the experimental results in Fig.2.The numerical results show a good agreement with the experimental one,thus the presented model for fracture propagation is verified.3.2.Optimal spacing for sequential fracture treatmentsAn example case is calculated in this section for spacing opti-mization in sequentially fractured horizontal well.The input pa-rameters used in this paper originated from a Chinese oilfield.The reservoir matrix is tight sandstone with low porosity and perme-ability.The pre-existing natural fractures are not considered in the simulations.The rock matrix is assumed to be isotropy and ho-mogeneous.The permeability and other formation parameters used in the simulations are the initial values of the reservoir matrix before fracturing.The parameters were selected approximately to keep business secret.The simulation of sequential fracturing for eight fractures was performed.The fracture spacing is set to be 30m,50m,70m and90m,respectively.The distance ofthe Fig.1.Schematic representation of fracture network andAFNW.parison offield measured surface pressure and the simulated one during fracturing.C.Liu et al./Journal of Natural Gas Science and Engineering29(2016)329e336 332exterior fractures to the model boundary is set to be400m to avoid boundary effects.The injection time is30min for a single fracture. The other parameters are listed in Table1.Due to the symmetry of geometry,material and load of the model,a half of the model is considered for all the simulations.As a fracture propagates,the in-situ stress contrast in the vicinity of the fracture is not a constant due to the variation of fracture net pressure(Wang et al.,2015b). The regions where in-situ stress contrast is less than3MPa during fracturing process is regarded as fracture network.The red(in the web version)contour shown in Fig.3indicates the in-situ stress contrast is less than3MPa.The fracture network produced by the first fracture is small due to the large initial in-situ stress contrast of the formation.Thefirst fracture propagates orthogonal to the horizontal well and its half-length is150m.The subsequent fractures grow away from thefirst one to some extent due to the stress perturbation induced by the propped fractures.The fracture length isfluctuating with each additional fracture.The fracture aperture is not a con-stant along fracture trajectory.The fracture width must be large enough to avoid sand plug when the proppant is injected into a fracture.The fracture widths at injection and deflection point of each fracture are regarded as an indicator for preventing sand plug in this section.The critical value of fracture width is2.55mm in practical engineering.At a small fracture spacing(30m,Fig.3a),the fracture width at the deflection point of the sixth fracture(D6)and the injection point of the seventh fracture(I7)is2.0mm and 1.1mm,respectively,which is less than2.55mm.So the small spacing(30m)is unacceptable for sequential fracturing in this case. The mechanical interaction between fractures decreases with fracture spacing.As the spacing varying from50m to90m,the fracture widths at injection and deflection points are larger than the critical value.Fig.4shows the AFNW with the varying fracture spacing.The AFNW increases with each sequential fracture for different spacing. The values of AFNW for the spacing of30m and50m do not show much difference.And the slops for the two curves are steep from the fourth fracture to the eighth.After the seventh fracture stim-ulated,the AFNW has nearly reached its maximum,and things do not change much for additional fractures.When the spacing is larger than50m,the AFNW decreases with the spacing.So extremely low or large spacing is not acceptable for horizontal well stimulation.Due to the risk of sand plug for the small spacing (30m),the optimal fracture spacing in this case is50m.3.3.Optimal spacing for simultaneous fracture treatmentsIn this section,the simulations of simultaneous fracturing for two and three fractures were presented.Optimal fracture spacing is achieved by calculating AFNW as well as fracture aperture with the varying fracture spacing.3.3.1.Optimal spacing for two fracturesThe propagation of two fractures with simultaneous injection method varying with four fracture spacing,20m,40m,60m and 80m are studied.The injection time is65min.The rock matrix permeability is0.1mD.The other simulation parameters are listed in Table1.The trajectory and fracture network for varying fracture spacing are shown in Fig.5.The fractures tend to grow away from each other due to the mechanical interaction between fractures. The fracture aperture along the propagation path at different fracture spacing is larger than the critical value(2.55mm),which satisfies the requirement for preventing sand plug.At the fracture spacing of20m and40m,the fracture network area between fractures is large,indicating the high mechanical interaction be-tween fractures.The AFNW with the varying fracture spacing is demonstrated in Fig.6.At a fracture spacing of20m,the AFNW is 25m,which is approximate to the spacing of40m(AFNW¼24m). In practical,the spacing of40m is preferable to20m due to the fewer perforation clusters in a specific horizontal well without significant production loss.The AFNW decreases dramatically with the fracture spacing after the spacing of40m is reached.At a spacing of80m,the AFNW is almost zero.So the optimal fracture spacing in this case is40m.3.3.2.Optimal spacing for three fracturesIn this study,the simulation of simultaneous fracturing for three fractures was performed.The fracture spacing is set to be30m,Table1Input parameters for simulation cases.Input parameters Value UnitsYoung's modulus50GPa Poisson's ratio0.2Fluid viscosity10mPa s Injection rate10m3/min Tensile strength2MPa Formation permeability0.5mD Fracture height80m Initial fracture half-length1m Initial pore pressure40MPa Maximum in-situ horizontal stress60MPa Minimum in-situ horizontal stress50MPa Porosity0.08Fig.3.Illustration of fracture aperture and fracture network at different spacing;the injection direction is from right to left.(a)Fracture spacing¼30m;(b)fracture spacing¼50m;(c)fracture spacing¼70m;(d)fracture spacing¼90m.Fig.4.Evolution of the AFNW with each sequential fracture for different fracture spacing.C.Liu et al./Journal of Natural Gas Science and Engineering29(2016)329e33633340m,50m,70m and 90m,respectively.Other parameters of this study are the same as previous two fractures propagation cases.The trajectory and fracture network for different fracture spacing are shown in Fig.7.Exterior fractures have greater lengths than interior fractures,and the length difference is due to the increasedcompression induced by fracturing.At a spacing of 30m,the frac-ture width of exterior fractures at the injection point is only 2.1mm,which is less than the critical value (2.55mm).So the spacing is unacceptable in the simultaneous fracturing process.When the spacing exceeds 30m,the fracture width is larger than the critical value due to the decreased interaction between fractures.At a spacing of 90m,the produced fracture network distributed only in the vicinity of the fractures with a small region.Fig.8shows the evolution of AFNW with the varying fracture spacing.The AFNW decreases dramatically with fracture spacing.However,the spacing of 30m is unacceptable in the fracturing process as discussed above.So the optimal fracture spacing in this case is 40m.3.4.Parameters effectsThis study analyzes the effect of Young's modulus,injection rate and in-situ stress contrast on fracture network and optimal spacing.The variations of the three parameters are not independentforFig. 5.Fracture network produced by two fractures with simultaneous injection method at different fracture spacing.(a)Fracture spacing ¼20m;(b)fracture spacing ¼40m;(c)fracture spacing ¼60m;(d)fracture spacing ¼80m.Fig.6.AFNW with the varying fracture spacing for two fractures with simultaneous injectionmethod.Fig.7.Fracture network produced by three fractures with simultaneous injection method at different spacing.(a)Fracture spacing ¼30m;(b)fracture spacing ¼50m;(c)fracture spacing ¼70m;(d)fracture spacing ¼90m.Fig.8.AFNW with the varying fracture spacing for three fractures with simultaneous injection method.C.Liu et al./Journal of Natural Gas Science and Engineering 29(2016)329e 336334。