Abstract

Two-phase (gas+water) flow is quite common in tight sandstone gas reservoirs during flowback and early-time production periods. However, many analytical models are restricted to single-phase flow problems and three-dimensional fracture characteristics are seldom considered. Numerical simulations are good choices for this problem, but it is time consuming in gridding and simulating. This paper presents a comprehensive hybrid model to characterize two-phase flow behaviour and predict the production performance of a fractured tight gas well with a three-dimensional discrete fracture. In this approach, the hydraulic fracture is discretized into several panels and the transient flow equation is solved by the finite difference method numerically. A three-dimensional volumetric source function and superposition principle are deployed to capture the flow behaviour in the reservoir analytically. The transient responses are obtained by coupling the flow in the reservoir and three-dimensional discrete fracture dynamically. The accuracy and practicability of the proposed model are validated by the numerical simulation result. The results indicate that the proposed model is highly efficient and precise in simulating the gas/water two-phase flow and evaluating the early-time production performance of fractured tight sandstone gas wells considering a three-dimensional discrete fracture. The results also show that the gas production rate will be overestimated without considering the two-phase flow in the hydraulic fracture. In addition, the influences of fracture permeability, fracture half-length, and matrix permeability on production performance are significant. The gas production rate will be higher with larger fracture permeability at the early production period, but the production curves will merge after fracturing fluid flows back. A larger fracture half-length and matrix permeability can enhance the gas production rate.

1. Introduction

Hydraulic fracturing technology can achieve the extraction of hydrocarbon from unconventional reservoirs at a low cost and improve the well performance significantly, which has been attracting considerable attention in recent years [17]. To generate highly conductive hydraulic fractures to increase well productivity, almost thousands of fracturing fluids are pumped into the reservoir [8]. The two-phase (gas+water) flow in the hydraulic fracture is quite common during the flowback and early-time production periods [912]. A common challenge for tight sandstone gas reservoir evaluation is the characterizations of the three-dimensional fracture. The productivity prediction of tight gas wells is greatly affected by the characteristics of the gas-water two-phase flow and fracture seepage parameters. Much research has been dedicated to theoretical models of fractured tight gas wells [1316], but most are only applicable to single-phase flow problems. Moreover, the three-dimensional fracture characteristics is seldom considered in these models.

Recently, flow behaviour analyses of flowback and early-time production from hydraulically fractured tight gas wells have been well known and researched. There is an array of analytical, semianalytical, and numerical simulation forecasting methods that can be used to analyze flowback data and make forecast for tight gas wells. Analytical models are based on the assumption that hydrocarbon does not breakthrough from fractured tight gas wells [1720]. However, these methods are limited to the single-phase flow of fracturing fluid in the hydraulic fracture and are generally inappropriate for analyzing flowback data after hydrocarbon breakthrough. Numerical simulation can rigorously consider complex flow physics, such as multiphase flow, reservoir, PVT properties, and fracture characteristics. It has gained much attention and been employed to investigate the two-phase flow behaviour and production performance of fractured tight gas wells during flowback [2124]. However, time-intensive computation performance and the difficulty of gridding make this method less attractive. In recent years, several semianalytical models have been developed to simulate two-phase flowback characteristics in tight oil and gas reservoirs [2527]. This approach can capture the two-phase flow of gas and water during the entire flowback period and have much higher computational efficiency than numerical simulation methods, which has attracted significant attention. In their model, the details of the transient two-phase flow in the fracture are incorporated. However, to simplify the mathematical model, the fluid flow in fractures is assumed as a 1D flow and the flux variation along the vertical direction in the fracture is ignored. Therefore, the semianalytical models cannot accurately capture the three-dimensional fracture characteristics, resulting in significant errors in production performance analysis.

This work aims at providing a hybrid model for simulating fracturing fluid flowback in tight sandstone gas wells considering a three-dimensional discrete fracture. We discretize the fracture into several panels, and each panel is treated as a plane source. The orthogonal transformation method and Green’s function method are employed to obtain the solution in the matrix system analytically. The finite-difference method is used to characterize the gas/water two-phase flow in the hydraulic fracture numerically. The solution of the model is semianalytically solved by coupling the transient flow in the fracture and matrix with the continuity of pressure and flux on the fracture surfaces. Based on the proposed model, the precision and efficiency of the proposed model are validated by using the commercial numerical simulator Eclipse. In addition, the influence of the two-phase flow in the hydraulic fracture on transient rate responses during flowback is analyzed. After that, the effects of several critical parameters on production performance are analyzed. The novelty of the proposed model is in the ability to semianalytically model the fracturing fluid flowback of fractured wells with three-dimensionally discrete fracture. The main advantage of the hybrid solution over more conventional numerical simulation is the reduced model setup and runtime without the loss of essential physics.

2. Methodology

2.1. Model Assumption

The scheme of a fractured tight sandstone gas well is shown in Figure 1(a). Figure 1(b) is the plan view of the well in the xy plane. The reservoir is box shaped, and the external boundaries of the top, bottom, and side are all assumed to be closed. The fracture is vertical, in the rectangle shape, and fully penetrate the formation. The fracture is initially filled with fracturing fluid after fracture stimulation and contributes to the fracture fluid recovery during flowback. The fluid flow in the fracture is considered as a 2D flow, and the gas enters the wellbore only through the fracture. This work assumes that water is only mobile in the fracture system to simplify the model. This physical problem is similar to what was assumed by Clarkson and Williams-Kovacs [11] and Jia et al. [28] to describe two-phase flowback behaviour in unconventional reservoirs. Some simplifying physical model assumptions for the derivation of the governing equation are listed as (1) the reservoir is assumed to be isotropic, (2) Darcy’s flow is considered both in the fracture and formation, and (3) the gas and water are produced with a constant bottom hole pressure.

2.2. Mathematical Model

As illustrated in the physical model, the hydraulic fracture system is the flow pathway directly connected with the wellbore. The fluids in the fracture system first flow into the wellbore, followed by the matrix-fracture flow. The fluid flow from the reservoir through the fracture and into the wellbore is composed of the two systems. This section will describe the mathematical model for the fluid flow in the hydraulic fracture system and in the matrix system.

2.2.1. Fluid Flow in the Hydraulic Fracture

Compared with the length and height of the fracture, the width of the fracture is minimal. So, the fluid flow in the three-dimensional fracture can be simplified as a two-dimensional flow and the mass conservation equations for the gas and water phases are given as
(1)ξρguξηρguη+q~~SCFg=tϕFρgSg,(2)ξρwuξηρwuη=tϕFρwSw,(3)Bc=ρSCρc,c=g,w
where uξ and uη are the darcy flow rates, ξ and η represent the direction of the fracture, Bg and Bw are gas and water volume coefficients, respectively, ρg and ρw are gas and water densities, respectively, Sg and Sw are gas and water saturations, respectively, ϕF is fracture porosity, and t is time.
Based on the Darcy law, the flow equation in the hydraulic fracture is written as
(4)uξ=βkFkrcμcpξ,(5)uη=βkFkrcμcpη,
where μg and μw are gas and water viscosities, respectively, krg and krw are relative permeabilities of the gas and water phases, respectively, p is the reservoir pressure, and β is the transmissibility conversion factor, where β=5.6146 for the field unit system.
The saturation relation equation is expressed by
(6)Sg+Sw=1.

Substituting equations (3) through (6) into (1) and (2), the governing equations in the three-dimensional fracture for the gas and water phases can be obtained.

For the gas phase, the governing equation can be described as
(7)ξβkFkrgSgμgpFBgpFpFξ+ηβkFkrgSgμgpFBgpFpFη+q~~SCFg+q~~SCWg=tϕFpFSgBgpF.
For the water phase,
(8)ξβkFkrwSgμwBwpFξ+ηβkFkrwSgμwBwpFη+q~~SCWw=tϕFpF1SgBw.
For the initial condition,
(9)pFξ,η,t=0=pi.
The well is assumed to produce at constant pressure, and the well production equation can be written as follows [29]:
(10)q~~SCWg=2πβkFkrgwFμgBglnreq/rWpFpW,q~~SCWw=2πβkFkrwwFμwBwlnreq/rWpFpW,(11)req=0.14ΔlF2+ΔhF2,
where kF is fracture permeability, pF is fracture pressure, pi is initial reservoir pressure, q~~SCFg is a source term representing the gas flux entering the fracture from the matrix at a point on the fracture surface per unit volume at standard conditions, q~~SCWg and q~~SCWw represent the gas and water production rate per unit volume at standard conditions, respectively, rW is the wellbore radius, wF is the width of the fracture, ΔlF is the length of the fracture panel, and ΔhF is the height of the fracture panel.

2.2.2. Fluid Flow in the Tight Sandstone Gas Reservoir

Based on the model assumption, for an isotropic-permeability reservoir, the governing differential equation for a single gas phase can be expressed as
(12)2pMx2+2pMy2+2pMz2+q~δxxw,yyw,zzw=ϕMμgpMctMpMkMpMt.
The initial condition is given by
(13)pMx,y,z,0=pi.
The external boundaries, including the top, bottom, and side, are assumed to be closed and given by
(14)pM0,y,z,tx=0,pMxe,y,z,tx=0,
for the side boundary,
(15)pMx,0,z,ty=0,pMx,ye,z,ty=0,
for the side boundary,
(16)pMx,y,ze,tz=0,
for the top boundary, and
(17)pMx,y,0,tz=0,
for the bottom boundary, where kM is the matrix permeability, pM is the matrix pressure, ϕM is the matrix porosity, ctM is the total compressibility of the matrix system, xe, ye, and ze are dimensions of the reservoir, (xw,yw,zw) is the center position of the plane source, and q~ is a point source influx from the matrix to the fracture.

2.3. Solution to the Mathematical Model

2.3.1. Solution of the Numerical Fracture Flow Model

In our work, the finite difference method is used to obtain the solution in the fracture system. The fracture is first discretized into small panels, as shown in Figure 2. Then, we apply the finite difference approximation to equations (7) and (8) at the time level n+1 and the implicit finite-difference form of fracture panel i,j can be obtained.

For the gas phase,
(18)lψi,jTgl,i,jn+1pFln+1pFi,jn+1+qSCFgi,jn+1+qSCWgi,jn+1=Cgpi,jn+1pFi,jn+1pFi,jn+Cggi,jn+1Sgi,jn+1Sgi,jn.
For the water phase,
(19)lψi,jTwl,i,jn+1pFln+1pFi,jn+1+qSCWwi,jn+1=Cwpi,jn+1pFi,jn+1pFi,jn+Cwgi,jn+1Sgi,jn+1Sgi,jn,
where i=1Nx and j=1Ny, NF is the total fracture panels and NF=NxNy, ψi,j is a 2D array whose elements are the existing neighboring segments of fracture panel i,j, and Tgl,i,jn+1 and Twl,i,jn+1 are the transmissibility of the gas and water phases between fracture panel i,j and its adjacent panels at time level n+1 and given by:
(20)Tgl,i,jn+1=Gl,i,jkrgl,i,jn+11μgBgl,i,jn+1,Twl,i,jn+1=Gl,i,jkrwl,i,jn+11μwBwl,i,jn+1.
The geometric transmissibility between fracture panel i,j and its adjacent panels is given by
(21)Gl,i,j=γlγi,jγl+γi,j,
where γi,j is the geometric transmissibility between fracture panel i,j and the interface, which is expressed as
(22)γi,j=βkFwFΔhFΔlFi,j.
The terms (Cgpi,jn+1, Cggi,jn+1, Cwpi,jn+1, and Cwgi,jn+1) in equations (18) and (19) are defined as
(23)Cgpi,jn+1=Vbi,jβΔtSgi,jnϕFi,jn+11Bgi,j+1Bgi,jnϕFi,j,Cggi,jn+1=Vbi,jβΔtϕFBgi,jn+1,Cwpi,jn+1=Vbi,jβΔt1Sgi,jnϕFi,jn1Bwi,j+1Bwi,jnϕFi,j,(24)Cwgi,jn+1=Vbi,jβΔtϕFBwi,jn+1,
where Vbi,j is the bulk volume of fracture panel i,j, Vbi,j=ΔlFΔhFwFi,j, the concepts of Copi,jn+1, Cwpi,jn+1, Cowi,jn+1, and Cwwi,jn+1 are similar to the compressibility coefficient, but are functions of pressure and saturation, and the coefficients of 1/Bgi,j, 1/Bwi,j, and ϕFi,j are given by
(25)1Bgi,j=1/Bgi,jn+11/Bgi,jnpFi,jn+1pFi,jn,(26)1Bwi,j=cwBw°,(27)ϕFi,j=ϕF°cF,
where Bw° is the water volume coefficient at a reference pressure, ϕF° is fracture porosity at a reference pressure, and cw and cF are the water compressibility and fracture compressibility, respectively.
Applying the finite difference approximation to equation (10) at the time level n+1, the implicit finite-difference form of production boundary condition can be expressed as
(28)qSCWgln+1=Tgl,Wn+1pFln+1pWn,lψW,(29)qSCWwln+1=Twl,Wn+1pFln+1pWn,lψW,
where pWn is the wellbore pressure.
The transmissibility of the gas and water phases between fracture panel l and the wellbore at the time level n+1 is given by
(30)Tgl,Wn+1=2πβkFwFlnreq/rWkrgln+11μgBgln+1,(31)Twl,Wn+1=2πβkFwFlnreq/rWkrwln+11μwBwln+1.
Combing equation (30) with (18), the gas flow equation can be rewritten as
(32)lψi,jTgl,i,jn+1pFln+1lψi,jTgl,i,jn+1+Cgpi,jn+1+Tgi,j,Wn+1pFi,jn+1+Cggi,jn+1Sgi,jn+1=qSCFgi,jn+1Tgi,j,Wn+1pWnCgpi,jn+1pFi,jnCggi,jn+1Sgi,jn
Substituting equation (31) into (19), the water flow equation can be rearranged as:
(33)lψi,jTwl,i,jn+1pFln+1+lψi,jTwl,i,jn+1+Cwpi,jn+1+Twi,j,Wn+1pFi,jn+1+Cwgi,jn+1Sgi,jn+1=Twi,j,Wn+1pWnCwpi,jn+1pFi,jnCwgi,jn+1Sgi,jn.
Equations (24) and (25) are the implicit finite-difference form of any fracture panels i,j at the time level n+1. Applying the two equations to NF fracture panels, the fracture flow equations can be arranged into a matrix format yield
(34)A×X=B+q,
where the matrix and vector in equation (34) are given by
(35)Ai,j=Tgi,j,Wn+1pWnCgpi,jn+1pFi,jnCggi,jn+1Sgi,jnTwi,j,WnpWnCwpi,jn+1pFi,jnCwgi,jn+1Sgi,jn,(36)X=pFi,jn+1Sgi,jn+1,Βl,i,jli,j=Tgl,i,jn+10Twl,i,jn+10,lψi,j,0000,else,Bi,j,i,j=lψi,jTgl,i,jn+1Cgpi,jn+1Tgi,j,Wn+1Cggi,jn+1lψi,jTwl,i,jn+1Cwpi,jn+1Twi,j,Wn+1Cwgi,jn+1,q=qSCFgi,jn+10.

2.3.2. Solution of the Reservoir Flow Model

Gas viscosity and gas compressibility factors in equation (12) are functions of formation pressure; to capture pressure-dependent properties and reduce the nonlinearity of the equations, a pseudopressure transformation is defined as follows
(37)m=20ppμgZdp.
And the pseudotime is defined as
(38)ta=0tμgictiμgp¯ctp¯dt,
where m is the pseudopressure and Z is the gas compressibility factor.
To make the equations homogeneous, some dimensionless variables are defined and tabulated in Table 1. Then, using the dimensionless variables and applying the pseudopressure transformation to equations (12)–(17) results in
(39)2mMDxD2+2mMDyD2+2mMDzD2+q~DδxDxwD,yDywD,zDzwD=mMDtD,(40)mMDxD,yD,zD,0=0,
for the initial condition,
(41)mMD0,yD,zD,tDxD=0,mMDxeD,yD,zD,tDxD=0,
for the side boundary,
(42)mMDxD,0,zD,tDyD=0,mMDxD,yeD,zD,tDyD=0,
for the side boundary,
(43)mMDxD,yD,zeD,tDzD=0,
for the top boundary, and
(44)mMDzD,yD,0,tDzD=0,
for the bottom boundary, where lR is the reference length, lR=ΔlF/2, mi is the initial pseudopressure of the reservoir, mM is the pseudopressure of the matrix system, and TM is the temperature.
In this section, the orthogonal transformation method is used to derive the analytical solution of the reservoir flow model. Applying the orthogonal transformation method to equations (39)–(44), one can obtain that
(45)λm~MD+q~DφxwD,ywD,zwD=m~MDtaD,m~MDtaD=0=0.
Equation (45) is the initial value problem of ordinary differential equations, and the solution of the equation can be expressed by
(46)m~MD=φxwD,ywD,zwD0taDq~DeλtaDτdτ,
where φ is a three-dimensional eigenfunction and λ is its eigenvalue.
The 3D eigenfunction and its corresponding eigenvalue can be obtained by solving the three 1D eigenfunctions and eigenvalues in the x, y, and z directions, respectively.
(47)2φxxD2=λxφx,φxxDxD=0=φxxDxD=xeD=0,(48)2φyyD2=λyφy,φyyDyD=0=φyyDyD=yeD=0,(49)2φzzD2=λzφz,φzzDzD=0=φzzDzD=1=0.
The one-dimensional eigenfunctions and eigenvalues of equations (47)–(49) are given by
(50)φxxD=cosαπxDxeD,λx=απxeD2,φyyD=cosβπyDyeD,λy=βπyeD2,φzzD=cosnπzD,λz=nπzeD2.
Thus, the three-dimensional eigenfunction and the corresponding eigenvalue can be expressed as
(51)φ=φxφyφz=cosαπxDxeDcosβπyDyeDcosnπzD,λ=λxλyλz=απxeD2βπyeD2nπzeD2.
After that, utilizing the inverse transformation of equation (46), the pressure solution for a three-dimensional volumetric source can be obtained and given by
(52)m~MDxD,yD,zD,tD=8q~DxeDyeD0taDα=01AxxcosαπxDxeDcosαπxwDxeDeαπ2τβ=01AyycosβπyDyeDcosβπywDyeDeβπ2τn=01AzzcosnπzDcosnπzwDenπ2τdτ.
Combing the definitions of the parameters in Table 1, we can rewrite equation (52) as
(53)mMmMx,y,z,t=αqSCTMϕMμgctMxeyeze0t1+2n=1expn2π2χtτxe2cosnπxwxecosnπxxe1+2n=1expn2π2χtτye2cosnπywyecosnπyye1+2n=1expn2π2χtτze2cosnπzwzecosnπzzedτ,
where mM is the current pseudopressure of the matrix system, mM is the initial pseudopressure of the matrix system, qsc is the withdrawal rate at standard condition, χ is the diffusivity and evaluated at the initial pressure, and α is the unit conversion factor.
The fracture is discretized into several fracture panels, and each fracture panel can be regarded as a plane source, as shown in Figure 3. Based on the theory of point-source function, integrating equation (53) along with the directions of the fracture panel results in
(54)mMmMx,y,z,t=αqSCTMϕMμgctMτ1τ2ΔlFi,j2ΔlFi,j2ΔhFi,j2ΔhFi,j21xe1+2n=1expn2π2χtτxe2cosnπxi,jxecosnπxxe1ye1+2n=1expn2π2χtτye2cosnπyi,jyecosnπyye1ze1+2n=1expn2π2χtτze2cosnπzi,jzecosnπzzedτdξdη,
where xi,j,yi,j,zi,j are the center position of the fracture panel i,j and equation (54) is the pressure response at an arbitrary position x,y,z at time t due to the contribution of a continuous plane source with a withdrawal rate of qsc.
According to the superposition principle, for the system in Figure 3, the pressure drop at the point (xo,yo,zo) caused by all fracture panels can be calculated and the equation is given by
(55)ΔmMxo,yo,zo,tn+1=mMmMxo,yo,zo,tn+1=αTMϕMμgctMk=1n+1i=1Nxj=1Nyq~SCMi,jktk1tkSxo,yo,zo,tn+1τ,xi,j,yi,j,zi,jdτ,
where q~SCMi,jk is the gas flow rate from the matrix to the fracture panel i,j per unit fracture area at time step k and S is the plane source function and expressed by
(56)Sxo,yo,zo,tn+1τ,xi,j,yi,j,zi,j=ΔlFi,j2ΔlFi,j2ΔhFi,j2ΔhFi,j21xe1+2n=1expn2π2χtτxe2cosnπxi,jxecosnπxxe1ye1+2n=1expn2π2χtτye2cosnπyi,jyecosnπyye1ze1+2n=1expn2π2χtτze2cosnπzi,jzecosnπzzedξdη.
To simplify the form of equation (55), let Ho,i,jn+1,k denote
(57)Ho,i,jn+1,k=αqSCTMϕMμgctMtk1tkSxo,yo,zo,tn+1τ,xi,j,yi,j,zi,jdτ.
Substituting equation (57) into equation (55) yields
(58)ΔmMxo,yo,zo,tn+1=mMmMxo,yo,zo,tn+1=k=1n+1i=1Nxj=1Nyq~SCMi,jkHo,i,jn+1,k.
With the superposition principle, the pressure response of fracture panel I,J caused by NF fracture panels at time level n+1 can be calculated by
(59)mMI,Jn+1=mMi=1Nxj=1Nyq~SCMi,jkHI,J,i,jn+1,n+1k=1ni=1Nxj=1Nyq~SCMi,jkHI,J,i,jn+1,k.
Equation (59) can be arranged into a matrix format and yields
(60)Hq~=m+mh.
In equation (60), the matrix H and vector h are defined as
(61)H=H1,1,1,1n+1,n+1H1,1,i,jn+1,n+1H1,1,NF,NFn+1,n+1HI,J,1,1n+1,n+1HI,J,i,jn+1,n+1HI,J,NF,NFn+1,n+1HNF,NF,1,1n+1,n+1HNF,NF,i,jn+1,n+1HNF,NF,NF,NFn+1,n+1,(62)h=k=1ni=1Nxj=1Nyq~SCMi,jkH1,1,i,jn+1,kk=1ni=1Nxj=1Nyq~SCMi,jkHI,J,i,jn+1,kk=1ni=1Nxj=1Nyq~SCMi,jkHNF,NF,i,jn+1,k.

And m=mM,mM,,mMT, m=mM1,1n+1,,mMi,jn+1,mMNF,NFn+1T, and q~=q~SCM1,1n+1,,q~SCMi,jn+1,q~SCMNF,NFn+1T.

2.3.3. Coupling Conditions for the Flow in the Two Systems

The flow in the matrix system and the discrete fracture system should be coupled to obtain the solution. Based on the flux exchange of the two systems and pressure continuity conditions on the fracture surface, the gas pressure continuity condition can be expressed as
(63)mFi,jn+1=mMi,jn+1.
And the gas flux continuity condition can be written as
(64)qSCFoi,jn+1=ΔlFi,jΔhFi,jwFi,jq~SCMi,jn+1.
We can group the fracture equations and the matrix equations to develop a system of equations to characterize the two-phase flow performance of the three-dimensional fracture.
(65)AX=B+q,Hq~=m+mh.

There are 3NF unknowns and 3NF independent equations in equation (65), and a unique solution is sure to be obtained. One should notice that the greatest challenge in solving equations (35) and (60) is to linearize equation (35). This is because the terms in matrix A and vector B are functions of pressure and saturation of the fracture panel, which are unknown parameters. This study uses the upstream weighting method for linearization in space and the simple iteration method for linearization in time. Then, the equations can be solved by a computer program. In the following, we analyzed the results based on this solution.

3. Results and Discussion

In this section, the proposed model is firstly validated using the numerical simulator Eclipse. Then, the effects of some critical fracture and formation parameters on rate transient behaviour during both the flowback period and the early-time production period are investigated.

3.1. Model Verification

To validate the model proposed in this paper, a fractured tight gas well model is built using the numerical simulator Eclipse and the schematic is shown in Figure 4. For comparison, the formation and fracture parameters are the same for both the Eclipse and the proposed model, which are listed in Table 2. The dimension of the whole reservoir is 600×800×10m, and the grid sizes for x, y, and z directions are 20 m, 20 m, and 10 m, respectively. The dimension of the fracture is 200×15×0.01m, and the grids near the hydraulic fracture face are gradually refined to capture the flow behaviour accurately. The fracture grids are initialized with two phases (gas+water), while the matrix is assigned 100% gas saturation. The gas PVT properties and the gas-water relative permeability curves in the hydraulic fracture system are shown in Figure 5.

The validation results are shown in Figure 6; it could be found that the result of the proposed model is in accordance with the numerical simulation result. The matching error is less than 2%. This shows the ability and precision of the proposed method in handling fracturing fluid flowback in tight sandstone gas wells considering three-dimensional discrete fracture. In addition, we also compared the computational speed of the two models. For the above simulations, the computation time for the numerical model is 48.2 s, while the computation time for the proposed model is 12.1 s. It is concluded that the proposed model is quite efficient in simulating the two-phase flow behaviour in three-dimensional discrete fractures.

3.2. Effect of the Two-Phase Flow in the Hydraulic Fracture

The main difference between the proposed model and the traditional models is that the three-dimensional fracture characteristics and two-phase flow behaviour are incorporated in the proposed model. In this section, we simultaneously simulated the rate transients of both the proposed and conventional models using the same group of formation and well parameters in Table 2. The initial water saturation in the fracture is given 0.9 and 0 in the proposed and conventional models, respectively. The modelling results are shown in Figure 7. It is shown that two-phase flow in the fracture system significantly influences the gas production rate at the early time. However, there is little influence after 60 days of production. In addition, the early gas production rate is smaller if the water-phase flow is considered in the fracture system. This is possible because the gas production from the tight matrix acts as a displacement effect for the limited water in the fracture system, so the two-phase flow in the fracture system only affects the gas production within the early flowback period. Besides, the gas production rate will be overestimated without considering the two-phase flow in the fracture.

3.3. Sensitivity Analysis

For tight sandstone gas reservoirs, the fracturing parameters (fracture permeability and fracture half-length) and even matrix permeability play a significant role in well performance and long-term production forecasting. Based on the proposed model in this paper, the influences of fracture permeability, fracture half-length, and matrix permeability on rate transient behaviour during the flowback period are analyzed. Except for the parameters analyzed, other parameters are the same and are shown in Table 2.

3.3.1. The Effect of the Fracture Permeability

The effect of the fracture permeability is analyzed in this part. Reference fracture permeability is given as 100, 300, 500, and 1000 mD, and the modelling results are shown in Figure 8. It should be noticed that gas and water exhibit higher rates with larger fracture permeability at the early production period. It is also observed that the fracture permeability significantly affects the gas rate during the early flowback period, while the curves will merge after 50 days of flowback. It is possible because the gas and water flow resistance in the fracture reduces with larger fracture permeability. In addition, the larger fracture permeability accelerates the water depletion in the fracture and the system will be dominated by the matrix earlier. So, the gas production rate becomes larger and the peak appears earlier as the fracture permeability increases.

3.3.2. The Effect of the Half-Length of the Hydraulic Fracture

Figure 9 illustrates the effect of the half-length of the hydraulic fracture on gas and water production rates. The half-length of the hydraulic fracture is given 80, 100, 120, and 150 m for the case study in this section. We could find that the effect of the length of the hydraulic fracture on the gas and water production rates is significant throughout most of the flowback period. As the length of hydraulic fracture increases, water exhibits higher production rates early but declines faster later. Besides, the gas production rate increases rapidly and the production peak appears later during the early flowback period with a larger half-length. This is possible because the half-length of hydraulic fracture reflects the relative capacity of fracturing fluid stored in the fracture system; a bigger half-length is the response of relative abundant reserves in the fracture system. For the fracture with a large length, the fracturing fluid at the far end generates strong support for water saturation near the wellbore. That is the reason why the gas production peak appears later as the fracture half-length increases.

3.3.3. The Effect of the Matrix Permeability

Figure 10 shows the effect of the matrix permeability on gas and water production rates. Herein, the matrix permeability is given 0.01, 0.05, and 0.1 mD for the case study and noted that matrix permeability significantly affects the gas production rate. The gas production rate increases significantly with a larger matrix permeability, and the production peak appears earlier at the early time. However, the matrix permeability has little effect on the water production rate during the flowback period. This is possible because the water flow behaviour is determined by the properties of the fracture system, while the gas flux from the matrix to the fracture and the matrix permeability dominate the associated gas production rate. Because the mobile water is assumed only in the fracture system, so, the effect of the matrix permeability on water production is not significant. In addition, the gas flow resistance in the matrix reduces with larger matrix permeability, which generates a high gas rate and strong support for fracture pressure. Consequently, the water production rate increases after the gas saturation develops in the fracture.

4. Conclusions

This work introduces and successfully applies a hybrid model to characterize two-phase flowback and early-time production behaviour in tight sandstone gas wells considering a three-dimensional discrete fracture. Using the hybrid approach, we investigate the influences of fracture permeability, fracture half-length, and matrix permeability on rate transient behaviour. The main conclusions of this work are as follows:

  • (i)

    Comparable reservoir simulator results reveal that our proposed model can be used to model gas and water flowback behaviour and early-time production performance of the three-dimensional fracture with satisfying precision and high computation efficiency

  • (ii)

    The fracture permeability, fracture half-length, and matrix permeability can significantly influence the rate transient behaviour. The gas production rate will be much larger with higher fracture permeability in the early flowback period, while the curves will merge after 50 days of flowback. The fracture half-length and the matrix permeability play a significant role in the gas production rate during the entire flowback and early-time production periods. A larger fracture half-length and matrix permeability can enhance the gas production rate

  • (iii)

    The gas production from the tight matrix acts as a displacement effect for the limited water in the fracture system, so the two-phase flow in the fracture system only affects the gas production within the early flowback period. However, the gas production rate will be overestimated by ignoring the two-phase flow in the hydraulic fracture

Nomenclature

     
  • Bg:

    Gas formation volume factor, Sm3/m3

  •  
  • Bw:

    Water formation volume factor, Sm3/m3

  •  
  • cF:

    Fracture compressibility, MPa−1

  •  
  • ctM:

    Total compressibility of the matrix system, MPa−1

  •  
  • cw:

    Water compressibility, MPa−1

  •  
  • Gl,i,j:

    Geometric transmissibility of fracture panel l and panel i,j

  •  
  • h:

    Formation thickness, m

  •  
  • kF:

    Permeability of the fracture, mD

  •  
  • kM:

    Permeability of the matrix, mD

  •  
  • krg:

    Relative permeability to gas, dimensionless

  •  
  • krw:

    Relative permeability to water, dimensionless

  •  
  • lR:

    Reference length, m

  •  
  • NF:

    Number of panels of the discretized fracture

  •  
  • pF:

    Fracture pressure, MPa

  •  
  • pi:

    Initial formation pressure, MPa

  •  
  • pM:

    Matrix pressure, MPa

  •  
  • pW:

    Wellbore pressure, MPa

  •  
  • qSC:

    Withdrawal rate at standard condition, m3/d

  •  
  • q~~SCFg:

    Gas flow rate entering the fracture from the matrix per unit volume at standard conditions, 1/d

  •  
  • q~~SCWg:

    Gas production rate per unit volume at standard conditions, 1/d

  •  
  • q~~SCWw:

    Water production rate per unit volume at standard conditions, 1/d

  •  
  • req:

    Radial distance of the side external boundary, m

  •  
  • rw:

    Wellbore radial, m

  •  
  • Sg:

    Gas saturation, dimensionless

  •  
  • Sw:

    Water saturation, dimensionless

  •  
  • t:

    Time, day

  •  
  • TM:

    Formation temperature, K

  •  
  • Tgl,i,j:

    Transmissibility of the gas phase between fracture panel l and i,j

  •  
  • Twl,i,j:

    Transmissibility of the water phase between fracture panel l and i,j

  •  
  • Tgl,W:

    Transmissibility of the gas phase between fracture panel l and the wellbore

  •  
  • Twl,W:

    Transmissibility of the water phase between fracture panel l and the wellbore

  •  
  • Vbi,j:

    Bulk volume of fracture panel i,j, m3

  •  
  • x,y,z:

    x, y, and z coordinates

  •  
  • xe,ye,ze:

    Reservoir dimension, m

  •  
  • xw,yw,zw:

    Center position of the plane source

  •  
  • ΔhFi,j:

    Height of fracture panel i,j

  •  
  • ΔlFi,j:

    Length of fracture panel i,j

  •  
  • Δt:

    Timestep size, d.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors acknowledge that this study was partially funded by the CNOOC Research Institute Co. Ltd. Science and Technology Projects (no. CNOOC-YXKJ-QZJC-TJ-2020-01).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).