Multi-fracture interactions during two-phase flow of oil and water in deformable tight sandstone oil reservoirs

Yongjun Yu, Wancheng Zhu, Lianchong Li, Chenhui Wei, Baoxu Yan, Shuai Li

Center for Rock Instability and Seismicity Research, Northeastern University, Shenyang,110819, China

ABSTRACT Tight oil reservoirs are complex geological materials composed of solid matrix,pore structure,and mixed multiple phases of fluids, particularly for oil reservoirs suffering from high content of in situ pressurized water found in China. In this regard, a coupled model considering two-phase flow of oil and water, as well as deformation and damage evolution of porous media, is proposed and validated using associated results, including the oil depletion process, analytical solution of stress shadow effect, and physical experiments on multi-fracture interactions and fracture propagation in unsaturated seepage fields. Then,the proposed model is used to study the behavior of multi-fracture interactions in an unsaturated reservoir in presence of water and oil. The results show that conspicuous interactions exist among multiple induced fractures.Interaction behavior varies from extracted geological profiles of the reservoir due to in situ stress anisotropy. The differential pressures of water and that of oil in different regions of reservoir affect interactions and trajectories of multi-fractures to a considerable degree. The absolute value of reservoir average pressure is a dominant factor affecting fracture interactions and in favor of enhancing fracture network complexity.In addition,difference of reservoir average pressures in different regions of reservoir would promote the fracturing effectiveness. Factors affecting fracture interactions and reservoir treatment effectiveness are quantitatively estimated through stimulated reservoir area.This study confirms the significance of incorporating the two-phase flow process in analyses of multifracture interactions and fracture trajectory predictions during tight sandstone oil reservoir developments.

Keywords:Multi-fracture interactions Two-phase flow Porous media deformation Hydraulic fracturing Continuum damage mechanics

1. Introduction

Hydraulic fracturing has been widely used as an efficient technique for reservoir stimulation, aiming at promoting oil and gas exploitation, particularly in shale, tight sandstone oil reservoirs,and other unconventional formations (Bunger et al., 2013a; Guo et al., 2015a; Meng et al., 2018; Zhao et al., 2019a). For the purpose of increasing production and stimulated reservoir volume(SRV), multistage clusters with perforations are preset along the wellbore prior to hydraulic fracturing. Multiple induced fractures emanate from the perforations and propagate into the target stratum.

Competition among multi-fractures is inevitable and influences the trajectories and propagation lengths among multi-fractures,and the interaction phenomena (stress shadow effect) have yet to be understood (Peirce and Bunger, 2015). In view of negative effects, interactions of multi-fractures constrain the initiation and propagation of some fractures during hydraulic fracturing. The fracture aperture is reduced by interactions, particularly for fractures in the interior of one cluster,leading to increasing perforation friction and entry resistance of subsequent injected fluid flowing into the fracture and further decreasing proppant conductivity(Weng,1993; Bunger et al., 2013b; Vahab and Khalili, 2018; Zhou et al., 2018; Cheng and Bunger, 2019). Meanwhile, the unexpected non-uniform flow rates partitioned among perforations are formed(Zhao et al., 2016; Cheng and Bunger, 2019). In contrast, the diverted fractures contribute to connecting oil and gas bearing areas that are originally unscheduled, thus enhancing the range of resource development and SRV, from the perspective of positive effects (Olson and Dahi-Taleghani, 2009; Kresse et al., 2013; Guo et al.,2015b).Fracture trajectory is an intuitive index for observing multi-fracture interactions.

Many factors, such as the reservoir pressure, intervals of perforations, clusters, and wellbores, injection rate, in situ stress anisotropy, and fracturing technique (e.g. simultaneous, sequential, or alternative fracturing), affect fracture trajectory and mechanical behavior of hydraulic fracture interference. The interference manifests diverse forms, such as mutual repulsive,sub-parallel and gradually coalescent growth of multi-fractures.Meanwhile, fracture interaction makes the fracture trajectory complex and affects the stress state of the rock mass (Pan et al.,2014; Zhang et al., 2014; Miao et al., 2018). Rabaa (1989) conducted physical experiments on fracture interference using hydrostones, and found that fracture geometry was influenced by the azimuth between horizontal wellbore axis and minimum in situ stress direction. Also, the repulsive propagation of multifractures was detected. Crosby et al. (2002) conducted experiments on two hydraulic fracture interactions, in which two independent injection tubes were used for individual pumping. The breakdown pressure of the second fracture was close to that of a previously induced fracture, indicating that stress shadow can be alleviated by appropriate arrangement of wellbores in certain stress regimes. A pioneer laboratory study on induced fracture propagation in unsaturated fields has been conducted by Bruno and Nakagawa (1991), in which the non-uniform seepage field was implemented by previously injecting at one perforation for a period of time, prior to formal fracturing. The fracture path was revealed to be severely affected by the previously injected perforation, due to its influence on effective stress alteration and strain energy release rate.

The tight sandstone oil reservoir should be regarded as a compound structure, composed of rock matrix, pores, and pressurized water and oil.Thus,two-phase flow process by co-seepage of water and oil has been considered to produce a complete description of reservoir performance (Pan et al., 2013; Asadi and Ataie-Ashtiani,2015; Guo et al., 2015b; Rozhko, 2016). Two-phase flow is a sophisticated process because the flow regimes are determined by a range of factors,including relative permeability of each phase fluid,and capillary pressure characteristics of the rock matrix. Examination of the influence of water relative permeability on production has found that increasing the water relative permeability favors higher gas production rates and water backflow (Cao et al., 2017).For a gas-water based two-phase flow system,simulation has been performed on coproduction of coalbed methane and shale gas in compound reservoirs (Meng et al., 2018). Also, the intersection mechanism of hydraulic and natural fractures has been discussed with extended finite element method (XFEM) results (Luo et al.,2018). However, these two simulations neglected capillary pressure, which reflects tortuosity of pore throat and other microproperties of pore structures (Van Genuchten,1980).

A system of partial differential equations pertinent to multiphase flow is always nonlinear and robustness-unsteady.Designing computational algorithms that are acceptable and versatility-capable has been of great interest (Asadi and Ataie-Ashtiani, 2015; Yan et al., 2018). Depending on the alternative choices of primary unknowns, two-phase flow coupled equations can be derived in distinct forms,which are evaluated based on the maximum iteration numbers and calculation execution time(Asadi and Ataie-Ashtiani, 2015). Pressure terms are regularly chosen as the primary unknowns to be implicitly solved rather than other configurations (e.g. pure saturation terms, or a mixed form in terms of pressure and saturation), where convergence is steadily approached in most cases with lowest solution oscillations (Asadi and Ataie-Ashtiani, 2015). Based on this, several hybrid models have been used to simulate the two-phase flow problems. Solid mechanics were solved by finite element method(FEM), XFEM, finite difference method and discontinuous displacement method (DDM), while flow regimes were determined by finite volume method(Ran et al.,1997;Lecampion et al.,2017).

Numerical simulations on multi-fracture interactions have been performed using various methods, such as XFEM, DDM and continuum damage mechanics (CDM) (e.g. Valkó and Economides,1994; Zhu and Tang, 2004; Bunger et al., 2012; Wu et al., 2012;Zhu et al., 2013a, b, 2016; Shojaei et al., 2014; Li et al., 2016, 2017;Li and Zhang, 2017; Luo et al., 2018). Luo et al. (2018) used the XFEM to study interactions among hydraulic and natural factures while considering two-phase flow of water and gas,and found that natural fractures were reactivated in different forms. Subsurface non-isothermal conditions have been considered during hydraulic fracturing, in which multiple branched fractures were induced by block crushing due to sudden thermal variations (Li et al., 2016,2017; Li and Zhang, 2017). Wu et al. (2012) investigated fracture interactions by unconventional fracture model(an enhanced DDM)in naturally fissured formations, and found that complex induced fractures resulted from the guidance effect by random natural fracture networks. CDM modeling based on FEM has some advantages, e.g.efficiently reproducing complete propagation process of hydraulic fractures (Valkó and Economides, 1994; Zhu and Tang,2004; Zhu et al., 2013a, b, 2016; Shojaei et al., 2014). In addition,accurate determination of stress singularity at the crack tip and tedious mesh refinement can also be avoided. In early times, the CDM based hydraulic fracturing simulations oversimplified the evolution of fracture length and width, where determination of locally zigzag fracture trajectories and dynamic changes of pore pressure of a heterogeneous reservoir was not available(Mu et al.,2019,2020;Valkó and Economides,1994;Shojaei et al.,2014;Zhao et al., 2019b).

The aforementioned numerical models concerned multiphase flow processes after reservoir stimulation treatments(Odeh,1981;Asadi and Ataie-Ashtiani, 2015; Guo et al., 2015b; Rozhko, 2016;Cao et al.,2017;Meng et al.,2018;Yan et al.,2018).To the best of the author’s knowledge, rare simulation or modeling was reported on the fracture propagation during a two-phase flow process. In addition, physical experiments implementing this complicated coupled process have remained unfeasible.

In this study, interactions of multiple hydraulic fractures are investigated in consideration of fracture propagation in unsaturated seepage fields. A fully coupled model incorporating twophase flow of oil and water as well as porous media deformation and damage evolution is proposed. The relative permeability and capillary pressure of two-phase fluid are also included. Elaborated validations are performed to demonstrate the versatility and capability of the proposed model. Afterwards, parametric studies are carried out to discuss the influence of reservoir average pressure changed by each phase fluid pressure on the multifracture interaction behavior and fracture propagation trajectories. In the end, the hydraulic fracturing effectiveness is quantitatively estimated based on stimulated reservoir area (SRA). The coupled model is based on CDM, in which SRA is determined equivalently by damage zone area(Zhu and Tang,2004;Zhu et al.,2013a, b, 2016).

2. Governing equations and verifications

Assumptions are made herein prior to modeling: (1) Hydraulic fracturing is an isothermal process such that thermal stress is neglected;(2)The plane strain assumption is adopted so as to allow discussion of two-dimensional (2D) hydraulic fracturing; (3) Tensile stress and strain are regarded as positive throughout this study;and(4)The reservoir pores are occupied by two kinds of fluids,e.g.water and oil, neglecting any dissolved gas phases in oil reservoir.

2.1. Coupled fluid-stress-damage (FSD) model considering twophase flow of oil and water

A set of coupled equations considering reservoir deformation and two-phase flow of oil and water is established. Damage evolution is coupled with porous media deformation (Zhu and Tang,2004; Zhu et al., 2013a, b). To consider reservoir heterogeneity,the Weibull function is used to generate mechanical parameters(e.g. Young’s modulus and uniaxial compressive strength),expressed as (Zhu and Tang, 2004; Zhu and Wei, 2010; Zhu et al.,2013a, b, 2016):

where χ denotes the arbitrary parameters,χ0is the average level of χ, andmis the homogeneity index describing the shape of distribution function and revealing the data fluctuation in vicinity of χ0.A largermrepresents a more homogeneous oil reservoir.

2.1.1. Mechanical equilibrium equation

In the coupled model, all the pertinent variables are averaged over the representative element volume (REV), as shown in Fig.1.

The elastic equilibrium equation of oil reservoir is expressed as

where σij(i,j=x,y,z) is the stress tensor andFi(i=x,y,z) is the body force.The strain-displacement relation is expressed as

where εij(i,j=x,y,z)is the strain tensor anduioruj(i,j=x,y,z)is the displacement vector. Reservoir average pressurepavcan be estimated as the averaged pressure of oil and water,defined as(Ran et al.,1997; Ma et al., 2017a; Luo et al., 2018; Yan et al., 2018):

Fig. 1. Schematic for the REV adopted in the coupled FSD model, where REV is decoupled into four subsections,including reservoir rock matrix,inner pore structures,and pressurized water and oil. The superposition principle is available for determination of REV’s total stress regimes.

Fig. 2. Schematic of the elastic damage evolution in uniaxial loading conditions.

whereSwandSoare the water and oil saturations,respectively;andpwandpoare the water and oil pressures, respectively. The water pressurepwis superimposed by the original water pressure in the reservoir, with the seepage pressure by water subsequently injected. Based on Eq. (4), the effective stress tensoris written as

where δijis the Kronecker’s delta equal to 1(wheni=j)or 0(wheni≠j);and α is the Biot’s coefficient,expressed as α = 1-Kb/Ks,in whichKbandKsare the bulk moduli of the porous media and the solid skeleton, respectively. Based on the isothermal deformation assumption, the constitutive relation for poroelastic media is expressed as (Zhou et al.,1998):

whereGis the shear modulus,ν is the Poisson’s ratio,and εVis the volumetric strain. Combining Eqs. (2), (3), (5) and (6), the twophase flow coupled equilibrium equation can be rewritten as a Navier’s equation form:

2.1.2. Damage evolution equations

The elastic damage evolution is adopted to simulate the gradual propagation of hydraulic fractures(Fig.2)(Zhu and Tang,2004;Zhu et al., 2013a, b, 2016).

The maximum tensile stress and Mohr-Coulomb criterion are used to decide the tensile and shear failures,respectively,described as

whereF1andF2are the two damage threshold functions denoting tension and shear failures, respectively; σ1and σ3are the maximum and minimum principal stresses,respectively;ft0andfc0are the uniaxial tensile and compressive strengths, respectively;and θ is the internal friction angle.The Young’s modulusEdecreases monotonically during the damage evolution process and is expressed by

whereE0is the initial Young’s modulus with no damage; andDis the damage variable ranging from 0 (intact state) to 1 (absolutely damaged state), and determined by(Zhu et al., 2013a; b, 2016):

where ε1and ε3are the maximum and minimum principal strains,respectively;εt0and εc0are the strains corresponding toft0andfc0,respectively;andnis a damage-based constitutive coefficient and is assumed to be 2 in this study (Zhu et al., 2013a; b, 2016). Notably,tensile failure is always judged predominantly over the shear failure in any loading scenarios. To consider damage effect on compactness of porous media,the Biot’s coefficient is rewritten as a step function:

where α0is the initial Biot’s coefficient. When the REV is completely damaged, regardless of tensile or shear damage mode,the Biot’s coefficient would reach a maximum value (i.e.1).

2.1.3. Two-phaseflow coupled equations

The masses of water and oil are described as

wheremwandmoare the water and oil masses, respectively; φ is the reservoir matrix porosity; and ρwand ρoare the water and oil densities, respectively. Based on mass balance equation of a single phase fluid (Zhu et al., 2016; Wang et al., 2017), mass balance equation of an oil-water system is expressed as

whereqwandqoare the flow rate vectors of water and oil,respectively;andQswandQsoare the source or sink terms(positive for source) of water and oil, respectively. Combining Eqs. (12) and(13) yields the fluid continuity equations for water and oil,respectively, and expressed as

The fluid compression coefficients of water and oil are defined as

whereCwandCoare the compression coefficients of water and oil,respectively; andVwandVoare the volumes of water and oil,respectively. The densities of water and oil change with pressure,which is described by(Bear,1972):

where ρwiand ρoiare the reference densities corresponding to reference pressurespwiandpoi,respectively.In this study,we refer to the standard atmospheric pressure 0.101 MPa where related ρwiand ρoidata are used.Based on Eqs.(16a)and(16b)and,time rate of density change is expressed as

The motion of the rock matrix particles can be depicted from the Euler’s viewpoint, where the mass conservation of rock is expressed as

where ρsis the rock density,which can be considered as a constant;andVsis the velocity vector of a rock mass point. Thus, Eq. (18) is further simplified as

Eq.(19)reflects the relationship between rock particle absolute velocity and matrix porosity. Based on the transport theorem of fluid dynamics,the substantial derivative of an arbitrary scalar(·)is expresses as (Bear,1972):

where D/Dtis the substantial derivative operation, andVis the velocity vector.Combining Eqs.(19)and(20),the following relation is met:

When considering the small deformation assumption, and neglecting inertia effect, the term accounting for migration accelerationV·?(·) in Eq. (20) can be neglected, which leads to simplification of the substantial derivative to the partial derivative as D(·)/Dt≈?(·)/?t, with acceptable accuracy. Hence, Eq. (21)becomes

The divergence ofVsin Eq. (22) is further expanded as

whereUis the displacement vector(U=ui+vj+wk);u,vandware the displacement components alongx-,y-andz-direction multiplied with base vectorsi,jandk,respectively.Substituting Eq.(23) into Eq. (22) yields

The reservoir matrix strain affects the porosity evolution during the hydraulic fracturing process, and thus we consider the change of porosity as (Tan et al., 2014):

where φ0is the initial porosity and it implies that the porosity reduction/increase results from matrix overall shrinkage/expansion. Substituting Eq. (25) into Eq. (24)yields

The constitutive equations describing the two-phase fluids are expressed by the modified Darcy’s law as

wherekmis the intrinsic permeability of the reservoir matrix;krwandkroare the relative permeabilities of water and oil,respectively;andgis the gravitational acceleration vector. In addition, the evolution of porosity and damage has an effect onkm,which is described by the cubic law multiplied by a damage-based exponential function,expressed as(Zhu and Wei,2010;Zhu et al.,2013a,b):

wherekm0is the initial intrinsic permeability of reservoir matrix;αkis the damage-permeability coefficient and is empirically set to 5(Zhu and Wei, 2010; Zhu et al., 2013a, b); andDand φ can be determined by Eqs. (10) and (25), respectively. Combining Eqs.((14),(17),(26)and(27),the governing equation of two-phase flow of oil and water, coupled with reservoir deformation and damage during hydraulic fracturing, can be expressed as

During derivation of Eqs. (29a) and (29b), the potential divergence operation, ?·[arb(x，y，t)g] = ?[arb(x，y，t)g]/?y, should be used and understood to account for gravitational effect only in vertical direction(i.e.y-axis in this study),in which arb(x，y，t)is an arbitrarily chosen field variable. As four unknowns are involved in the coupled two-phase flow equations (i.e.Sw,So,pwandpo), two other supplementary equations are needed for a closed-form solution of the proposed model.

The relative permeability is the key parameter in describing the relative seepage characteristics of multiple phases of fluid flowing through the porous media. Determining the relative permeability expression through a pure analytical solution is not feasible;instead, some empirical relations between relative permeability and wetting saturation were used(Pinder and Gray,2008;Cao et al.,2017; Ma et al., 2017a). In this study, we refer to an empirical exponential model illustrated in Pinder and Gray (2008), together with the experimental results of relative permeability curves of tight sandstone with lithology similar to the reservoir in this study.The sandstone samples are subjected to a two-phase seepage process of water and oil(Ren et al.,2018).The relative permeability model is expressed as (Pinder and Gray, 2008):

whereSwiis the irreducible water saturation;kmaxpis the oil phase saturation atSwi; andmandmpare two laboratory parameters to be determined, respectively. The fitting results show that parametersmandmpare equal to 0.85 and 4.67, respectively(Fig. 3).

Fig. 3. Relative permeability curves of water and oil measured by Ren et al. (2018), as shown in black lines. Curve fitting is carried out based on relative permeability model(Eq. (30a) and (30b)) in this study.

Fig. 4. Capillary pressure curve of sandstone measured by the MICP method (Zhang et al., 2018).

Fig.5. Oil-water capillary pressure curve converted from Zhang et al.(2018),and data are fitted to the VG model to provide an explicit relationship between capillary pressure and water saturation for performing numerical simulations.

The summation of water and oil saturations should be 1 at each time point, which gives the auxiliary equation for saturation as

When the interfacial tension between water and oil phases is incorporated,the capillary pressurepcrepresenting the differential pressure between the two-phase fluids is expressed as

A conventional capillary pressure model proposed by van Genuchten (VG model) is used to link the wetting saturation with capillary pressure with three pending coefficients to be fitted,expressed as (Van Genuchten,1980)

wherehcis the hydraulic head of the capillary pressure,defined ashc=pc/(ρwg);α0,NandMare the fitting parameters; andSweis the effective saturation of wetting phase, defined as

Combining Eqs. (33) and (34), we have

whereNandMare the experimental parameters to be decided.

Mercury injection capillary pressure(MICP)method is regularly used for capillary pressure measurement.For example,Zhang et al.(2018)carried out an experiment to characterize the pore structure of the oil sandstone by MICP method, and the capillary pressure curve is shown in Fig. 4.

For field applications,conversions are made when determining relationships of capillary pressure curves among different twophase flow systems (Rozhko, 2016; Yang, 2016). All data on the capillary pressure curve obtained by MICP method should be converted to an oil-water based two-phase flow system based on(Yang, 2016):

Fig. 6. Evolution of porosity during oil depletion with initial porosity φ0(x，y，0) = 0.25 (?x， y ∈AΩ).

Fig. 7. Comparison of the results from present model and Ran et al. (1997).

where subscripts “R” and “L” denote the reservoir and laboratory conditions, respectively; and σ and θ are the surface tension and contact angle,respectively.The surface tension of mercury is σHg=480 (m N)/m and the contact angle is θHg= 140°, while the surface tension of oil-water is σow= 25 (m N)/m, and the contact angle is θow= 0°(Yang, 2016). Based on Eq. (36), the capillary pressure of oil-water system (pow) in reservoir condition has relationship with capillary pressure (pHg) by MICP method aspow≈pHg/15, and the converted curve of capillary pressure of the oilwater system is shown in Fig. 5. Meanwhile, the VG model is used for curve fitting, with the obtained coefficients α0,NandMequal to 1.64 ×104, -0.46 and -12.77, respectively.

As the two-phase flow controlled equations are highly nonlinear, the physical time domain needs to be small enough,aiming at mitigating the oscillation process during nonlinear solutions.

2.1.4. Boundary and initial conditions

The boundary conditions of stress and displacement are prescribed for the mechanical solution as

Fig.8. Schematic of the experiment and simulation on hydraulic fracturing in unsaturated seepage fields in Bruno and Nakagawa(1991):(a)Sketch of the configured sample,and(b) FEM model established in this study.

Table 1Primary input parameters for fracturing simulations with non-uniform pore pressure fields.

(1) The Dirichlet condition:

(2) The Neumann condition:

Fig.9. Comparison of experimental and numerical results in scenario I.For the evolution of damage zone(b),we use“-1”denoting complete tensile failure,“1”denoting complete shear failure,and “0”denoting undamaged state.Values among the three values denote a transitional state.Same explanations are for other figures of damage zone evolution in this study.

Fig.10. Comparison of the experimental and numerical results in scenario II.

Fig.11. Schematic of a pressurized fracture in an infinite plate. Blue color represents the pressurized fluid providing a uniform distributed pressure P in the fracture, and a small red square represents a volume body close to the fracture.

whereu0(x，0), σ0(x，0) andp0(x，0) are the initial displacement,stress and fluid pressure within the solution domain Ω,respectively.

As the processes of porous media deformation and two-phase flow are fully coupled, the equations above are complicated and nonlinear,leading to no feasible analytical solutions.The proposed model,together with prescribed boundary and initial conditions,is implemented into Comsol Multiphysics, and the hydraulic fracturing process is carried out by further MATLAB compiling. With regard to two-phase flow equations, pressures are chosen as the primary unknowns and solved implicitly, in favor of a stable convergence and less solution oscillations, as discussed in Introduction.

2.2. Model validations

Elaborated model validations are performed covering numerical solution, analytical solution, and physical experiments. It includes verification of full branches of the coupled model by determining porosity evolution during oil exploitation in a water-containing reservoir (Odeh, 1981; Ran et al., 1997). Also, portions of the full branches when the flow equations are degenerated back to singlephase flow of water or oil are compared with stress shadow analytical solutions, laboratory results of fracture propagation in unsaturated fields(Bruno and Nakagawa,1991),and multi-fracture interactions during sequential fracturing (Bunger et al., 2011).

2.2.1. Validation against the porosity evolution during oil depletion

It is recommended that the model validation should be conducted firstly at a scenario where all the branches of the proposed coupled equations can be thoroughly verified. The porosity evolution during an oil exploitation process has been investigated in Ran et al. (1997), considering the coupled effects of reservoir deformation with two-phase flow of oil and water, with geological and flow parameters in Odeh (1981). Input parameters are as follows:Reservoir size of 3048 m×3048 m with a wellbore of 0.0762 m in radius; Young’s modulus of 1300 MPa; Poisson’s ratio of 0.15;flowing bottom hole pressure of 5 MPa; initial porosity of 0.25;permeability of 100 mD; initial reservoir pressure of 30.3 MPa;initial oil saturation of 0.88;and initial water saturation of 0.12.As the fracture propagation process is not considered in this validation, the damage variableDcan be regarded as zero to reflect a coupled process of solid deformation and two-phase flow,referring to the definition ofDin Eq.(10).For the related parameters above,notably that all engineering units adopted in Ran et al. (1997) and Odeh(1981)have been transformed as above.The porosity φ(x，y，t)is averaged over the entire domain to produce an overall porosityevolution as

whereAΩis the area of the pay zone.The overall porosity evolution during the oil depletion process from the production wellbore is obtained using the proposed model (Fig. 6).

Oil depletion causes the reservoir matrix shrinkage, thus reducing overall reservoir porosity with production time. A comparison of the simulation result with that of Ran et al.(1997)for the overall reservoir porosity is shown in Fig. 7.

Both the present model and Ran et al. (1997)’s model predict a similar descending trend in porosity evolution with limited distinctions. Therefore, the proposed model is competent in reproducing the fully coupled process of reservoir matrix deformation with the two-phase flow of oil and water.

Fig.12. Distribution of Young’s modulus(Pa)with different homogeneity indices(m).A central line connecting upper surface of fracture and top edge of plate is selected to record induced stresses, as denoted by blue dotted line AB. Point A is in the middle of upper surface of fracture, and point B is in the middle of top edge.

Fig.13. Comparison of induced stresses obtained by the analytical solution (Sneddon and Elliot,1946) and the proposed model.

Fig.14. Schematic for (a) the experiment on multi-fracture interactions and (b) the corresponding FEM model for validation.

2.2.2. Validation against fracture propagation in unsaturated pore pressurefields

Sandstone slabs with three perforations were used as the rock samples to investigate the influence of preexisting non-uniform pore pressure fields on the propagation trajectory of the fractures induced thereafter (Bruno and Nakagawa,1991). The schematic of the experiment is illustrated in Fig.8(Bruno and Nakagawa,1991).

Table 2Brief input parameters for validation against experimental results on fracture interactions.

A square rock slab ofL=15.2 cm is well constrained with steel bolts. Three perforations of 5 mm in diameter are numbered withA1,A2andA3, respectively. The notch is 0.6 mm in aperture and 6 mm in length. In scenario I, perforationA2is pre-injected with a kind of mineral oil to establish a constant fluid pressurePa2= 1.73 MPa prior to fracturing. This formed non-uniform pore pressure is maintained for 5 h. Then, the upper perforationA1is injected with the same fluid at a rate ofQa1=2 cm3/min,providing an increasing fluid pressurePa1over time. Throughout the experiment, the fluid pressurePa3in perforationA3is maintained at ambient pressure, lower thanPa1andPa2. In scenario II, both perforationsA1andA3have the same pre-injected pressure of 2.07 MPa, which forms another non-uniform pore pressure distribution pattern.For model validation against these experiments,the conditionsSw=krw=pw=0 andSo=kro=1 are met.In addition,the body force termFiin Eq. (7) should be regarded as ρog, accounting for fluid weight during the vertical seepage process.Input parameters are briefly listed in Table 1.

The fracture initiation and propagation in scenarios I and II are shown in Figs.9 and 10,respectively,with the experimental results for comparisons (Bruno and Nakagawa,1991).

Fig.15. Numerical results of multi-fracture interactions in comparison with experimental results (Bunger et al., 2011) at the end of pumping: (a) Multi-fracture interactions and propagation trajectories; (b) Comparison between experimental and numerical results at the end of pumping. The legend shows the Young’s modulus in Pa.

Regardless of the distinct patterns of preexisting pore pressure before fracturing, the propagation of induced fractures in both numerical simulations and experiments show the similar trajectories, demonstrating the dominant role of pore pressure fields in diverting fracture propagation to the region with higher pore pressure. From these validations (Figs. 9 and 10), we confirm the significance to investigate multi-fracture interaction behavior in unsaturated seepage fields,such as a two-phase flow process of oil and water in the actual oil reservoir status.

2.2.3. Validation against the stress shadow effect

During hydraulic fracturing, the stress regime in vicinity of the induced fracture would change because of the additional stresses provided by the pressurized fluid.The extra stresses are commonly denoted as“stress shadow”effects.A landmark study has proposed an analytical solution of induced stresses around a pressurized Griffith crack in an infinite plate of free constraint (Sneddon and Elliot, 1946); the schematic of this analytical model is shown in Fig.11.

The analytical solution of induced normal stresses σii,inducedis expressed by(Sneddon and Elliot,1946):

Fig.16. Curve of pumping pressure versus injection time in numerical simulation.

wherePis the uniform distributed pressure of injected fluid, and other geometrical relations are shown in Fig.11.

In comparison to the analytical solution, a rectangular plate of 1 m×1 m provided with a center-penetrated fracture is modeled.The fracture is 0.1 m in length and 0.005 m in aperture.The length ratio of fracture to the plate edge is 1/10, to eliminate boundary effects.Nonetheless this is only an approximation to the analytical infinite model. The proposed flow equations are simplified such that only water flow is considered. Thus, the conditionsSw=krw=1 andSo=kro=po=0 are met.Input parameters include Young’s modulus (10 GPa), Poisson’s ratio (0.3), uniform water pressureP(10 MPa), water compressibility (5 × 10-10Pa-1),permeability (1 × 10-30mD), Biot’s coefficient (1 × 10-30), water viscosity (1 mPa s), and water density (1000 kg/m3). Both permeability and Biot’s coefficient are set to be infinitesimal,such that the time-dependent seepage pressurePof fluid is neglected and only the mechanical effect is considered, in accordance with the analytical solution (Sneddon and Elliot, 1946). Heterogeneous plates are modeled by settingm(see Eq. (1)) to be 100 (nearly homogeneous) and 3 (heterogeneous), respectively. The Young’s modulus distribution is shown in Fig.12.

Fig. 18. Model established for Case Study 1 (plan view): (a) Schematic of the model with geometrical sizes and perforation configurations,where scales are distorted for a clear identification; and (b) Sketch of FEM model with the perspective position of a wellbore and load conditions.

Fig.17. Multi-fracture interactions when the sample size is enlarged 2-fold to 700 mm×700 mm compared with that used in Bunger et al.(2011).Yellow dashed square outlines the original scope in the experiment (Bunger et al., 2011) as shown in Fig.15.

Fig.19. Model established for Case Study 2 (front view): (a) Schematic of the model with sizes and lithology, where geometrical scales are distorted for a clear identification; and (b) Sketch of FEM model with the positions of wellbores and load conditions.

A comparison of the numerical and analytical solutions for induced normal stresses σηη,induced and σξξ,induced along lineABis shown in Fig.13.

Stress results predicted by numerical simulation were found to coincide well with the analytical solutions (Fig. 13a and b). More fluctuations occur whenm=3.With the distance to middle pointAincreases, the induced stresses decrease dramatically to zero,reflecting that the induced stresses (or stress shadow effect) have limited influential areas in the intermediate vicinity of a pressurized fracture.

Bunger et al.(2011)conducted laboratory tests on interactions of multiple hydraulic fractures by means of sequential fracturing of four stages.The cubic gabbro samplewas 350 mm×350 mm×350 mm in size and was characterized as a kind of granite.The schematic of theirexperiment and numerical model established in this study are shown in Fig.14.The spacing between each perforation was 0.015 m,and the four staged fracturing had a sequence of 1 → 2→3 → 4(from bottom to top).The fracturing fluid was a water solution combined with blue dye and glycerin. In this validation, the mixed injected solution is regarded as wetting phase withSw=krw=1 andSo=kro=po=0;other input parameters are listed in Table 2.

Table 3Primary input parameters for Case studies 1 and 2.

Computational steps 1—100, 101—200, 201—300 and 301—400 are set to represent fracturing stages 1—4, respectively, such that the steps 100, 200, 300 and 400 denote the end time of each fracturing stage.The results are shown in Fig.15,with 1 step equal to 5 s.

Longer and shorter fractures are observed to emerge alternately from bottom to top perforations during sequential fracturing,with the propagation paths in experiment and simulation coincident.Four primary fractures numbered 1—4 are formed around the four slotted notches. Fractures 1 and 3 propagate mainly along the direction of horizontal stress σx;while fractures 2 and 4 deviate from the expected direction,due to the interference by fractures 1 and 3.Model validation again confirms that the preexisting pressurized fractures impose extra stresses on the adjacent fractures, which suffer from diversion from the preferable propagation direction.The curve of pumping pressure versus injection time of simulation is shown in Fig.16.

The breakdown pressures of the four stages are 21.2 MPa,21.4 MPa,22.2 MPa and 21.4 MPa,respectively,compared with the experimental results of 25.1 MPa,23.3 MPa,25.9 MPa and 26.9 MPa,respectively (Bunger et al., 2011). Lower breakdown pressures are obtained in the numerical simulation,which is probably due to the heterogeneous differences between the numerical model and real rock sample. It has been verified experimentally and numerically that a shallow fracture near the rock surface tends to be in a saucer shape, thus reflecting boundary effects on the path of fractures(Bunger et al., 2013a). Therefore, it makes sense to discuss boundary effects by increasing the sample size in Fig.15. In this respect,another numerical simulation is conducted in which the model size is magnified 2-fold compared to that in Fig.15, with other conditions and parameters keeping unchanged. The numerical results are shown in Fig.17.

Fig. 20. Numerical results of the Base case in Case Study 1.

Overall fracture interactions are seen to still exist when the sample size is enlarged, in which fractures 1 and 3 could propagate farther than fractures 2 and 4 do. The interaction behavior is similar to that shown in Figs. 14 and 15. However,when provided a larger space for fracture propagation, fracture 1 still propagates to the edges of the model, which in turn retards the propagation of fractures 2, 3 and 4. Thus, boundary might have affected some local fracture propagation, while maintaining approximately unchanged the behavior of multi-fracture interaction.

Fig. 21. Gradual change of the direction of maximum compressive effective principal stress negative) during injection process for the Base case in Case Study 1. Red arrow denotes direction and the length of each arrow is proportional to the local size of

Table 4Different initial water and oil pressures in Cases 1 and 2 of Case Study 1(unit:MPa).

3. Numerical simulations

Models for simulating multi-fracture interactions in tight oil reservoir are established while considering the two-phase flow of oil and water.Subsurface oil reservoir is a combination of different rock strata,where the sandstone acts as a pay zone,and the stratum mainly composed of mudstone acts as enclosing barriers sandwiching the pay zone. Case Study 1 simulates the situation of a horizontal geological profile of the pay zone, coplanar with the horizontal wellbore axis. Case Study 2 simulates the situation of a longitudinal geological profile containing pay zone and barriers.The longitudinal geological profile is perpendicular to the horizontal wellbore axis,thus considering a sandwiched reservoir.The input physical parameters for both Case studies 1 and 2 are collected from the log data of the Shengli Oil Field in Shandong Province, Eastern China, which belong to typical oil blocks in a clearly water-containing state. The reservoir has characteristics of low permeability, complex tortuosity of pore structures, and abundant reservoir water capacity, making it a typical unconventional reservoir of tight sandstone.

3.1. Model setup for Case Study 1

A pay zone of 1 m×1 m is established in a horizontal geological section,with three perforations slotted on a perspective horizontal wellbore.Schematic of the model and the meshed FEM model with boundary configurations are shown in Fig.18.

Fig. 22. Numerical results of Case 1 in Case Study 1.

Fig. 23. Numerical results of Case 2 in Case Study 1.

Initial water and oil pressures of the pay zone are denoted bypw0_pandpo0_p, and prescribed aspw0_p= 0.4 MPa andpo0_p=1.6 MPa,respectively.The sum of water and oil pressures is 2 MPa.Thus,reservoir average pressurepavis 1.42 MPa,calculated by Eq.(4).This configuration of initial conditions is regarded as the“Base case” in Case Study 1.

3.2. Model setup for Case Study 2

Barrier rock layers composed of mudstone have relatively higher strength and lower porosity and permeability than the pay zone,which is very common in the Shengli Oil Field. Tight argillaceous rock contentis in awide range of 6.1%—63.9%,derived fromfield logs.Barriers and pay zone have different distributions and contexts of water and oil.As can be seen in Fig.19,two perforated wellbores are separated horizontally with a distance of 20 cm. Three points denoted withA,BandCare set,to record the variations of effective stress components in orthogonal direction,i.e.PointAis located at the center of pay zone,and pointsBandCare located at the centers of the upper and bottom barriers,respectively.

Initial water pressures of the pay zone and barriers are denoted bypw0_pandpw0_b, and initial oil pressures in the pay zone and barriers are denoted bypo0_pandpo0_b, respectively. The initial pressures are set aspw0_p= 0.4 MPa,po0_p= 1.6 MPa andpw0_b= 1.6 MPa,po0_b=0.4 MPa.From this,the sum of water and oil pressures is 2 MPa both for pay zone and barriers; while the reservoir average pressures of pay zone and barriers are 1.42 MPa and 0.95 MPa, respectively, calculated using Eq. (4). This configuration of initial conditions is regarded as the “Base case” in Case Study 2. The input parameters for Case studies 1 and 2 are briefly listed in Table 3.

For both Case studies 1 and 2, all perforations are fractured simultaneously by a monotonic water pressure increment Δpw= 0.75 MPa per step, where 1 step is equal to 10 s of seepage calculation.The coefficient of in situ stress differenceKσis used to estimate in situ stress anisotropy as follows:

Table 5Different initial water and oil pressures in Cases 3—6 of Case Study 1 (unit: MPa).

Fig. 24. Initial distributions of water and oil pressures in Cases 3—6 of Case Study 1. All legends are in unit of Pa.

Fig. 25. Fracture interactions in Cases 3—6 of Case Study 1. The legend shows the Young’s modulus in Pa.

(1) Case Study 1:

(2) Case Study 2:

As for Case Study 2,σ3should be regarded as the average of in situ stresses applied on the right boundary of reservoir (Eq. (42b)and Fig.19). It shows that in situ stress anisotropy in Case Study 2 is much lower than that in Case Study 1.

4. Results and discussion

4.1. Multi-fracture interactions in Case Study 1

4.1.1. Base case

Fig. 20 shows the simulation results of the Base case. The fractures propagate mainly along the maximum in situ stress direction without any apparent deviations. The fracture emanated from the left perforation propagates with a higher rate, compared with fractures emanated from the middle and right perforations,due to local heterogeneity difference among areas in the vicinity of the three perforations. The middle fracture is the shortest since it is subjected to extra compressive stresses (stress shadow effect)provided by both the left and right fractures. Injected water forms an elliptical conjunctive area around the three perforations, with relatively high water saturation (Fig. 20d) and low oil saturation(Fig. 20e), reflecting a process in which injected water drives the preexisting oil in the pay zone. Mutual interference among the hydraulic fractures is due to local in situ stress alterations,and the gradual diversion of the maximum compressive effective principal stressis determined (Fig. 21).

Fig. 26. Fracture trajectories when only middle perforation (a) or right perforation (b) exists. The legend shows the Young’s modulus in Pa.

4.1.2. Influence of water and oil pressures within the pay zone

The influence of different allocations of water and oil pressures on multi-fracture interactions is investigated in two cases, where the sum of water and oil pressures is kept as 2 MPa(Table 4).

The reservoir average pressurepavis determined by Eq. (4),which provides a difference ofpavbetween the north and south parts of the pay zone. For Cases 1 and 2 of Case Study 1,pavis close to that of the Base case.The water and oil have the opposite trends of convection based on the provided initial two-phase fluid pressures, while the direction of potential convection is just aligned with preferred fracture propagation direction (i.e.yaxis). The results of Cases 1 and 2 are shown in Figs. 22 and 23,respectively.

The two exterior perforations can be fractured more effectively to achieve longer induced fractures than that of the Base case(Fig.20).Meanwhile,the degree of fracture propagation diversion is increased,with apparent repelling propagation observed in the two exterior fractures, and the middle perforation is still suppressed from effective fracturing.The existing pressure difference of water phase and that of oil phase enhance convective flow of the two immiscible fluids. Each phase fluid in the pay zone regions with higher pressure has the potential to flow to other pay zone regions with lower pressure,such that the overall enhancement of the twophase convective flow is formed in Cases 1 and 2 of Case Study 1.Once the perforations are effectively fractured with induced fractures, the injected fluid seeps into the pay zone and then immediately participates in the two-phase convective flow,which in turn promotes overall seepage motion in the pay zone and increases effective tensile stress till failure.The time consumed for hydraulic fracturing is less for Cases 2 and 3 of Case Study 1 than that of the Base case, when comparing computational step with fracture length.

Another four cases are prescribed,in which initial water and oil pressures are the same as the Base case in some regions of the pay zone, while other regions have a declined initial pressure for both water and oil,as listed in Table 5.Reservoir average pressurepavis determined by Eq. (4). These initial pressure configurations make both water and oil possess a trend of convection in the same direction, from west to east, while the direction of potential convection is orthogonal to the preferred fracture propagation direction. In Cases 3—6,pavis obviously distinct from that of the Base case,for a gradualpavdecrease from west to east regions of the pay zone. The initial distributions of water and oil pressures in Cases 3—6 are shown in Fig. 24.

The initial distributions of water and oil pressures in the pay zone can be clearly observed by combining Fig.24 with Table 5.The fracture interactions in Cases 3—6 of Case Study 1 are shown in Fig. 25.

Interactions of fractures are seen to be similar to those in the Base case(Fig.20),with predominance of the left fracture,while the middle perforation is the most constrained. The time required for fracturing is reduced from Case 3 to Case 6,compared with the Base case, according to the computational step with fracture length.Even though the convection direction is orthogonal to the preferred fracture propagation direction, the enhanced potential of relative seepage between water and oil also accelerates the fracturing propagation rate, compared with the Base case of uniform distributions of water and oil pressures.

4.1.3. Special case: A single perforation

In order to validate that both the middle and right perforations have the same ability to be fractured, two cases are selected in which only a single perforation exists(middle or right perforation),with other parameters being the same as that of the Base case.Fig. 26 shows the fracture trajectories when only a single perforation exists.

Fig. 27. SRA evolutions in all cases versus injection time, and initial absolute values of differential pav in Base case and Cases 1—6 of Case Study 1 referring to Tables 4 and 5.

Table 6Three SRA data in Case 2 for data anomaly analysis.

Fig. 28. Numerical results of the Base case in Case Study 2.

Fig. 29. Gradual change of the direction of maximum compressive effective principal stress negative) during injection process of the Base case in Case Study 2. Red arrow denotes direction and the length of each arrow is proportional to the local size of

Both the middle and right perforations are observed to be fractured effectively, forming an approximately equal length in each fracture,which verifies that all the three perforations have the same capability to produce effective induced fractures. Increasing SRV is the goal to maximize permeability enhancement, fracture complexity and ultimate production.Instead,SRA can be used in a 2D simulation (Ren et al., 2016). From the perspective of damage mechanics,the SRA is calculated from the damaged zone area,as an index for fracturing effectiveness evaluation after reservoir treatment.Fig.27 presents the changes of SRA in Base case and Cases 1—8 of Case Study 1, also the curves of SRA versus injection time normalized to that of the Base case.Normalized SRA reflects a clear variation of fracture propagation rate with respect to the Base case,based on which fracturing effectiveness can be quantitatively evaluated.

The SRA in Cases 1—6 of Case Study 1 is larger than that in the Base case(Fig.27a and c).More apparent distinctions are observed in Fig. 27b and d. Comparing the Base case with Cases 1 and 2(Fig. 27a, b and f), it suggests that the normalized SRA can be increased apparently by increasing the difference of reservoir average pressure in different regions of the pay zone. Comparing the Base case with Cases 3—6 (Fig.27c, d and f),it shows that the difference of reservoir average pressure in different regions of the pay zone can increase the normalized SRA,although the reservoir average pressure in Cases 3—6 is lower than that of the Base case(Table 5).The left perforation shows predominance in some cases(Figs. 20 and 25); however, Fig. 27e confirms the equivalent capability of hydraulic fracturing among three perforations. It should be pointed out that some abnormal fluctuations in Fig.27b and d are meaningless (in dotted ellipses), because they occur in positions prior to the actual ascending of SRA, as observed in Fig.27a—c.The values on curve of SRA versus injection time are all infinitesimal, before SRA increases apparently when induced fracture propagates unsteadily, and all these values should be regarded as zero.Taking Case 2 as an example,three SRA data are displayed in Table 6.

It can be observed that the abnormal data of normalized SRA stem from the first two rows of Table 6,which is likely to stem from calculation accuracy affected by iteration algorithms.Therefore,the abnormal large normalized SRA prior to the obvious increase of SRA should be neglected.

4.2. Multi-fracture interactions in Case Study 2

4.2.1. Base case

Fig. 28 presents the numerical results of the Base case in Case Study 2.The fracture curvature in Base case of Case Study 2 is larger than that of Case Study 1,by comparing Fig.28 with Fig.20,because the in situ stress anisotropy in Case Study 2 is lower than that in Case Study 1 (see Eq. (42)). Two induced fractures emanate from the perforations, and show no obvious interactions (Step 51). The changes of water and oil saturations are shown in Fig. 28d and e,which reflect increased water content and decreased oil content during hydraulic fracturing. With water being continuously pumped into the two wellbores, the fractures propagate fast in the combined rock layers and water has the propped effect on the fracture surfaces newly updated. By this time, the mutual interference among the hydraulic fractures becomes apparent,and leads to unstable fracture propagation.Interactions are seen as a kind of“repelling”behavior,making the initiated hydraulic fracture curves away from each other. Theoretically, the repelling effect can be explained by stress shadow effect according to Eq. (41). The pressurized hydraulic fractures transfer the fluid pressure from fracture surfaces into the surrounding rocks, thus enhancing the compressive in situ stress along the direction vertical to preferred propagation path and diverting the originally parallel fractures to propagate in the reoriented direction.

Fig. 29 shows the gradual change of the direction of maximum compressive effective principal stressAt the initial state,is mainly along longitudinal direction of maximum in situ stress σv(Step 1).Thengradually tends towards diverting into horizontal direction by continuous injection, especially in the pay zone,attributed to stress shadow effect of the induced fractures. The diversion becomes more obvious with further fracture propagation(Steps 55 and 59). The trend of maximum compressive effective principal stress mainly matches the fracture propagation (Fig. 28a and b and 29).

4.2.2. Influence of water and oil pressures between the pay zone and barriers

To understand how the fracture interactions are affected by different water and oil pressure distributions in the pay zone and barriers, several cases of initial water and oil pressures are prescribed (Table 7). The sum of water and oil pressures remains 2 MPa,the same as that of Base case.The reservoir average pressurepavis estimated by Eq.(4)and close to that of Base case.These initial configurations provide gradual increases ofandfrom Cases 1 to 3,and the water has the potential flow from barriers to the pay zone, while the oil has the potential flow from pay zone to barriers.

The fracture trajectories in Cases 1—3 of Case Study 2 are shown in Fig. 30. The general profiles of facture interactions remain the mutual repelling in Cases 1—3, which is similar to that of the Base case. Fracture propagation rates are improved a little from Case 1 to Case 3, by comparing fracture lengths with computational steps. The pressure difference of each phase fluid increases the relative seepage motion of that phase, thus enhancing the overall seepage motion of the immiscible two phase fluids in the sandwiched reservoir. The subsequent injectedwater into the wellbore not only fractures the rock, but also participates in the convective flow, leading to further enhancements of relative seepage motion.

Table 7Different initial water and oil pressures in Cases 1—3 of Case Study 2 (unit:MPa).

Another three cases are selected in which both water and oil pressures increase in the pay zone while maintaining the same two-phase fluid pressures in barriers as that of the Base case,producing a gradual increase of reservoir average pressure in the pay zone from Cases 4—6 of Case Study 2 (Table 8). The reservoir average pressurepavis determined by Eq.(4),and the difference ofpavbetween pay zone and barriers is more apparent in Cases 4—6 than that of the Base case.

Interactions of fractures in Cases 4—6 of Case Study 2 are presented in Fig. 31. Interactions of the fractures remain the mutual repelling as before, and the fracture propagation rate increases from Case 4 to Case 5, by comparing fracture lengths and steps.However,the propagation regime is really distinct in Case 6,where the mutual repelling tends to be the most apparent, leading to anearly horizontal diversion of the right upper fracture. Some complex transverse induced fractures emanate from the wellbores,and connect to make a larger fragmented region in the pay zone.Once the pay zone is fractured, a large amount of stored elastic energy is dissipated. The overall stiffness and strength of the pay zone decrease, leading to a greater difference in the mechanical properties between pay zone and barriers than those at the initially undamaged state. The penetrating propagation of the fractures from pay zone to barriers is restricted, due to the resistance from the barriers.Some fractures are diverted into propagating along the pay zone-barrier interface for some distance,before turning back to the preferable direction. Fracture diversion contributes to increasing induced fracture complexity, which is expected for larger SRV and oil-gas outputs. In general, the fracture diverted propagation along reservoir interface is realized by a lowerKσ(Eq.(42b)),a higher effective stress(Eq.(5))due to the increases of both oil and water pressures, as well as the preexisting of laterally configured rock interface, from the aspect of geometrical characteristic of the interbedded reservoir.

Table 8Different initial water and oil pressures in Cases 4—6 of Case Study 2 (unit:MPa).

Fig. 30. Interactions of induced fractures in Cases 1—3 of Case Study 2. The legend shows the Young’s modulus in Pa.

Fig. 31. Interactions of induced fractures in Cases 4—6 of Case Study 2. The legend shows the Young’s modulus in Pa.

Fig. 32. Changes of absolute value of differential effective in situ stresses at points A, B and C versus injection time.

The diverted propagation of fractures is caused by in situ stress reorientation.Based on Eqs.(4)and(41),in situ stress reorientation occurs in the region that satisfies the following condition:

Eq.(43)can be further rewritten in a differential effective in situ stress form:

Fig. 33. The fracture propagation trajectories with a single wellbore: (a) Only the left wellbore; and (b) Only the right wellbore. The legend shows the Young’s modulus in Pa.

The overall trend of in situ stress reversal effect of Case 6 is greater than that in the Base case,particularly for pointAlocated in the vicinity of two wellbores. Therefore, it confirms that the complex fracture network in Case 6 results from local reversal of effective in situ stresses,which can be regarded as a positive factor for volume fracturing.

4.2.3. Special case: A single wellbore

In order to validate that each fracture emanated from the two wellbores has the same ability to propagate,another two cases are designed to consider a single wellbore at the left or right of the reservoir profile, with other conditions being the same as that of the Base case. The fracture propagation trajectories are shown in Fig. 33.

In both cases 7 and 8, fractures propagate in the preferred direction of maximum in situ stress (σv), when there is no interference by another fracture. The right fracture approaches the same length as the left one,by 1 step(t=10 s)in advance,because of the local strength difference between the areas in the vicinity of the two wellbores. Fig. 34 shows the changes of SRA in the Base case and Cases 1—8 of Case Study 2, and also the curves of SRA versus injection time normalized to that of the Base case.

The SRA is noted to increase in Cases 1—6 of Case Study 2,compared with the Base case, as shown in Fig. 34a and c. More apparent distinctions are detected in Fig. 34b and d. Comparing the Base case with Cases 1—3 (Fig. 34a, b and f), it shows that even though the differences of reservoir average pressure between the pay zone and barriers are similar, the normalized SRA can also be increased, with increasing pressure difference of each phase fluid in different regions of the reservoir (referring to Table 7). Comparing the Base case with Cases 4—6 (Fig. 34c, d and f), the normalized SRA is increased to a larger magnitude than that of Cases 1—3, as both the reservoir average pressure and the difference of the reservoir average pressure are increased in Cases 4—6.Fig.34e again confirms the equivalent capability of hydraulic fracturing for left and right wellbores. The data fluctuation prior to obvious SRA increase can be explained as that in Case Study 1.

4.3. Discussion on competing propagation and two-phase flow parameters

Regarding to both Case Studies 1 and 2, several results reveal a predominant propagation of some fractures over others, e.g.Figs.20,25 and 28,where both model geometrical shape and load conditions are symmetrical. This is due to the local random distribution of heterogeneous mechanical properties (e.g. Young’s modulus and strength) of rock, which causes asynchronous initiation and propagation of fractures among different perforations.Once the asynchrony forms,it is impossible for the lagged fracture to pursue the previous one.The distinction of effectiveness among perforations was also found through physical experiment by Zhou et al. (2018), where two notches were slotted in the casing pipe to study the interference between two adjacent induced fractures in shale samples. It was also found in several cases that only one single obvious initiated fracture could be detected by means of acoustic emission and direct observation of the sample after pumping,reflecting the heavy restraint of one fracture on the other fractures during simultaneous fracturing.

Fig. 34. SRA evolutions in all cases versus injection time, and initial absolute values of differential pav in the Base case and Cases 1—6 of Case Study 2 referring to Tables 7 and 8

For a precise description of two-phase flow in fracture networks,the capillary pressure and relative permeability reflecting pore structural properties of the reservoir rock,i.e.the inner connection,tortuosity and diameter distribution of pore channels, and the relative seepage ability of two-phase fluids, need to be accurately considered (Ma et al., 2017b). However, in this study, the characteristics of relative permeability and capillary pressure between the pay zone and barriers are not distinguished, owing to the lack of available data for barriers, thus presenting a simplified simulation for the compound reservoir. The capillary pressure reflects the interfacial properties, and the interfacial tension determining the mechanical equilibrium between the two immiscible fluids of oil and water. The shape and estimation of the capillary pressure (pc-Sw) curve are varied among the reservoirs, resulting in that the solutions of fluid pressure are dependent on the capillary curve adopted. In addition, the relative permeability of oil and water is interrelated with the wetting phase saturation. Therefore, a chained relation ofkrw-kro-Sw-pcdoubtless has influence on the flow solutions.Furthermore,the pressure allocations for water and oil, and the overall reservoir pressure are prescribed as initial conditions for parametric studies. A detailed investigation lies in that fracture interactions and flow regimes are determined when the reservoir has been undergoing oil recovery for a period of time,where more complicated and diverse initial distributions of water and oil saturations, and capillary pressure should be provided as input initial conditions prior to numerical simulation. The above mentioned items, including the effects of mechanical property differences between the pay zone and barriers, as well as the rock interface mechanism and stress anisotropy, will be discussed in future, aiming at reflecting the indirect influences of reservoir characteristics and geological conditions on multi-fracture interactions during hydraulic fracturing.

5. Conclusions

In this study, the compound tight sandstone oil reservoir composed of solid matrix, pore structure, and pressurized water and oil is modeled to understand the interactions among multiple induced fractures during hydraulic fracturing.A set of FSD coupled equations considering two-phase flow of oil and water, porous media deformation and damage evolution processes is proposed.The coupled model is validated against several literature results,including the oil depletion process, analytical solution of induced stresses, and physical experiments on fracture interactions and propagation in unsaturated seepage fields, prior to numerical applications. Some main conclusions are derived as follows:

(1) Interactions among multiple induced fractures are conspicuous,both in vertical and horizontal geological profiles, and curved repelling behavior is more frequently observed in vertical profile of the sandwiched reservoir with lower in situ stress anisotropy than that in horizontal profile within the pay zone.The competing initiation and propagation of multifractures result in an asynchronous propagation process and different ultimate fracture regimes.

(2) The difference of reservoir average pressure examined by pressure difference of water, and that of oil in different regions of the reservoir are found to have an influence on the fracture trajectory and interaction behavior, and contribute to accelerating hydraulic fracturing. Promoted hydraulic fracturing effectiveness occurs even though the reservoir average pressures remain proximity.

(3) Fracture interactions and fracturing effectiveness are also sensitive to the absolute value of reservoir average pressure.The synchronous increased water and oil pressures in different regions of the reservoir promote the fracture mutual interference, as well as produce fracture network complexity. Differential pressures in different regions of reservoir, caused by increased absolute value of reservoir average pressure in some regions,play a more important role in raising the hydraulic fracturing effectiveness to a larger degree than those cases with similar reservoir average pressures.

(4) Fracturing effectiveness can be estimated intuitively through SRA for cases designed in this study,in which the fracturing effectiveness enhancements by increasing the difference of reservoir average pressure and the absolute value of reservoir average pressure are quantitatively analyzed. Fracture interaction behavior and fracturing effectiveness apparently suffer from local non-uniform distribution of water and oil pressures to a considerable degree, notably by referring to normalized SRA curves. Hence, it demonstrates the significance of considering the two-phase flow process for accurate analyses of fracture interaction behavior and fracture trajectory predictions.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interests associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

The present work is funded by National Natural Science Foundation of China (Grant Nos. 51761135102 and 51525402), and the Fundamental Research Funds for the Central Universities(Grant No.N180105029). These supports are gratefully acknowledged. The first author is grateful to Dr. Sanbai Li (Northeastern University,China)for the discussion and terminology calibration of this work.