# Production Decline Analysis for Two-Phase Flow in Multifractured Horizontal Well in Shale Gas Reservoirs

Authors: Wei-Yang Xie; Xiao-Ping Li; Lie-Hui Zhang; Xiao-Hua Tan; Jun-Chao Wang; Hai-Tao Wang
Journal: Journal of Chemistry
 (2015)

---

## Abstract

After multistage fracturing, the flowback of fracturing fluid will cause two-phase flow through hydraulic fractures in shale gas reservoirs. With the consideration of two-phase flow and desorbed gas transient diffusion in shale gas reservoirs, a two-phase transient flow model of multistage fractured horizontal well in shale gas reservoirs was created. Accurate solution to this flow model is obtained by the use of source function theory, Laplace transform, three-dimensional eigenvalue method, and orthogonal transformation. According to the model’s solution, the bilogarithmic type curves of the two-phase model are illustrated, and the production decline performance under the effects of hydraulic fractures and shale gas reservoir properties are discussed. The result obtained in this paper has important significance to understand pressure response characteristics and production decline law of two-phase flow in shale gas reservoirs. Moreover, it provides the theoretical basis for exploiting this reservoir efficiently.


---

## Body

## 1. Introduction

In order to obtain industrial gas flow, horizontal drilling and multistage hydraulic fracturing techniques have been widely used in SGR (shale gas reservoirs) development. Compared to the conventional reservoirs, SGR have their own uniqueness; for example, the shale is a self-sourcing reservoir in which the gas is stored either by compression or by adsorption on the surface of the solid material. Transport of gas in SGR is a complex multiscale flow, which is from macroscale to the molecular scale. Numerous researchers have studied migration mechanisms in SGR, but the two-phase flow caused by multistage hydraulic fracturing is usually ignored for analyzing and forecasting a shale gas well.According to Hill and Nelson [1], 20% to 85% of the total gas storage in SGR may be subject to sorption. Montgomery et al. [2] pointed out that the absorbed gas in SGR is mainly stored in pores, natural fractures, and matrix surface and desorption rate would change with the difference of desorption time, desorption pressure, and organic matter volume. Freeman et al. [3] and Cheng [4] employed a numerical simulator to study the flow regimes for a multistage hydraulic fractured horizontal well in SGR and the effect of gas desorption was briefly discussed. The studies of hydrofracture are mainly concentrated on the fractures’ characteristics and fluid flow mechanism. Xiao et al. [5] presented an analytical model for effective thermal conductivity of nanofluids by considering the effect of Brownian motion of nanoparticles. Tan et al. [6] presented a dual fractal reservoir transient flow model created by embedding a fracture system simulated by a tree-shaped fractal network into a matrix system simulated by fractal porous media and discussed the influence of different fractal factors on pressure transient responses. Cai et al. [7–9] analyzed the capillary-driven flow in gas-saturated porous media based on fractal theory and modified the classical Lucas-Washburn equation. Later, the hydraulic conductivity and effective porosity in porous media have been investigated by many researchers, such as Tan et al. [10]. Freeman et al. [11] presented a complicated numerical model to study the adsorption and desorption in SGR and described the change of gas composition during wells production. After that, some special fractal reservoir models have been presented by Tan et al. [12, 13]; in these papers they discussed dual fractal flow characteristic and liquid carrying phenomenon.Some scholars employed new methods to analyze production decline of SGR. Medeiros et al. [14] discussed the analysis of production data from hydraulically fractured horizontal wells in shale reservoirs. Ozkan et al. [15] presented a discussion of fractured horizontal-well performance in conventional and unconventional reservoirs and provided interpretation of the different production performance of fracturing horizontal wells in both types of formations. After that, scholars presented rate prediction models [16–18] in SGR, but the proposed models may have some limitations when applied to analyze production data in a shale well. Cho et al. [19] presented the results of an experimental study of pressure-dependent natural-fracture permeability in SGR and found that the effect of pressure-dependent natural-fracture permeability on shale gas well production was a function of the permeability of the matrix system. J. H. Lee and K. S. Lee [20] studied the multiscale heterogeneity in the existence of various permeability magnitudes and spatial distributions of permeability. Du et al. [21] presented a new experimental method explaining the flooding tests in the low-permeability reservoirs. Alkouh and Wattenbarger [22] presented a gas-water model to solve the challenge of simulating a natural fracture with trapped water without imbibitions in SGR. In the same year, Alkouh et al. [23] added post-frac water flowback and long-term water production to the analysis and showed the benefits of a new method for combining water flowback and long-term water production data in shale gas analysis. Adefidipe et al. [24] presented a new model to describe the two-phase flow in shale gas wells and analyzed the obvious gas phase coning phenomenon. Xie et al. [25] presented a transient flow model which considered two-phase flow after hydrofracture in shale gas reservoirs to analyze the transient pressure response.Until now, the studies of a multistage hydraulic fractured horizontal in SGR are mainly concentrated on single phase gas flow; the two-phase flow caused by flowback after hydrofracture in SGR does not attract much attention. In order to understand the effect of two-phase flow in SGR, a two-phase flow transient analysis model of multistage fractured horizontal well in SGR is presented, and the accurate solution to the flow model is obtained using source function theory [26, 27], Laplace transform, three-dimensional eigenvalue method, and orthogonal transformation. Type curves of pressure response and production decline are plotted and different flow regimes are identified. The influence of different factors on production decline performance is discussed.


## 2. Physical Model

Figure1 shows a horizontal well with hydraulic fractures in the SGR center. The physical model assumptions are as follows.(1)
The SGR are dual mechanism and with sealed boundary.(2)
Two-phase flow exists in fractures and horizontal well.(3)
All fractures are assumed to be equally spaced, to penetrate the formation completely, and to be perpendicular to the horizontal wellbore. Both the horizontal wellbore and the fractures are assumed to be infinitely conductive.(4)
Ignoring the impact of gravity and capillary forces.Figure 1
Two-phase flow in multistage hydraulic fractured horizontal well in SGR.


## 3. Mathematical Model

### 3.1. The Two-Phase Flow Governing Equations in Fracture System

Governing flow equations of two-phase flow in fractures under three-dimensional Cartesian coordinates can be written as follows.Gas flow equation in fractures is(1)∇ρgKKfrgμg∇pf-qm=∂ϕfSgρg∂t.Liquid flow equation in fractures is(2)∇ρlKKfrlμl∇pf=∂ϕfSlρl∂t.The transient diffusion rate of the shale gas from the shale matrixqm can be expressed as follows:(3)qm=G∂cm∂rm.The two-phase flow pseudopressure is employed and defined as follows:(4)mp=∫0pfρgKfrgμg+ρlKfrlμldpf.With (1), (2), and (4), the two-phase flow governing equation can be written as follows:(5)∇2m-qm=ϕfK∂Sgρg+Slρl∂t.Employing Taylor expansion to simplify (5), we can obtain the following equation:(6)∇2m-qm=η∂m∂t,where (7)η=ϕfCtKρgKfrg/μg+ρlKfrl/μl.Employing source function into (6) yields(8)∇2m+1K-qm+qgpscRTsc+qlρl∇2m+1KqgpscRTsc·δx-xi,y-yi,z-zi,t-τ=η∂m∂t.To simplify the mathematical model, dimensionless parameters are defined as follows.The dimensionless two-phase pseudopressure is(9)mD=Khqtmi-m.The dimensionless wellbore storage coefficient is(10)CD=CKhL2ηρgKfrg/μg+ρlKfrl/μl.The dimensionless time is(11)tD=CDtL2η.The dimensionless coordinate value ofx is(12)xD=xL.The dimensionless coordinate value ofy is(13)yD=yL.The dimensionless coordinate value ofz is(14)zD=zL.The dimensionless relative mass concentration is(15)cmD=1ρscci-cm.The dimensionless radial distance in matrix system is(16)rmD=rmR.Using the dimensionless parameters, (8) can be written as follows:(17)∇2mD-hL2qtGρscR∂cmD∂rmD+qgpscRTsc+qlρl∇2mD-hL2qtGρscR·δxD-xiD,yD-yiD,zD-ziD,tD-τD=CDη∂mD∂tD.The two-phase flow mathematical model in Laplace space is obtained by taking the Laplace transformation [28] of (17). The two-phase flow mathematical model is as follows:(18)∇2mD¯-hL2qtqgpscRTsc+qlρlδxD-xiD,yD-yiD,zD-ziD=sCDηmD¯+hL2GρscRqt∂cmD¯∂rmD.With the condition of sealed boundary, the dimensionless form of the initial condition is(19)lim|tD|→0mD¯=0.The dimensionless forms of the outer boundary conditions are(20)lim|xD|→0∂mD¯∂xD=0,lim|xD|→a/L∂mD¯∂xD=0,lim|yD|→0∂mD¯∂yD=0,lim|yD|→b/L∂mD¯∂yD=0,lim|zD|→0∂mD¯∂zD=0,lim|zD|→h/L∂mD¯∂zD=0.


### 3.2. The Transient Diffusion Rate Equations in Matrix System

With Fick’s second law and the assumption of spherical matrix blocks, the transient diffusion rate equation can be written as follows [29]:(21)1rm2∂∂rmDrm2∂cm∂rm=∂cm∂t.Using dimensionless parameters and Laplace transform to rewrite (21),(22)∂2cmD¯∂rmD2+2rmD∂cmD¯∂rmD=R2uCDDL2ηcmD¯.For transient diffusion case, the center of the spherical matrix block can be treated as a no-flow boundary, so the inner boundary condition can be written as follows:(23)∂cmD∂rmDrmD=0=0.The concentration on the external surface of matrix block is in equilibrium with the pressure in the fracture. The outer boundary condition can be written as follows:(24)cmDrmD,tDrmD=1=cmDm.The general solution of (22) can be expressed as follows:(25)cmD¯=ArmDeWrmD+BrmDe-WrmD.With the boundary condition and general solution expression, the accurate solution of the dimensionless mathematical model of matrix system for transient diffusion case can be obtained as follows:(26)cmD¯=cmDm¯2shR2sCD/DL2ηrmDcmD¯=·eR2sCD/DL2ηrmD-e-R2sCD/DL2ηrmD.With the definition of hyperbolic functions, (26) can be rewritten as(27)cmD¯=cmDm¯shR2sCD/DL2ηrmDshR2sCDDL2ηrmD,where cmD(m)¯ is the dimensionless equilibrium volumetric gas concentration. Employing Langmuir isotherm [30], there is(28)cmDm¯=-VLmLmi-mfmL+mfmL+mi¯cmDm¯=-qtKhVLmLmL+mfmL+mimD¯.Equation (28) can be rewritten as(29)cmDm¯=-qtKhVLmLmL+mfmL+mimD¯=-αmD¯,(30)α=qtKhVLmLmL+mfmL+mi.Inserting (29) into (27) and getting derivative with respect to rmD, the solution of the matrix model in the Laplace domain is as follows:(31)∂cmD¯∂rmD=-αmD¯R2sCDDL2ηcoth⁡R2sCDDL2η-1.


### 3.3. The Two-Phase Flow Mathematical Model for Multistage Fractured Horizontal Well in SGR

Inserting (31) into (18) results in(32)∇2mD¯-hL2qtqgpscRTsc+qlρlδ=sCDη-hL2GρscRqtαR2sCDDL2ηcoth⁡⁡R2sCDDL2η-1·mD¯.To simplify, (32) can be rewritten as(33)∇2mD¯-hL2qtqgpscRTsc+qlρlδ=fsmD¯,(34)fs=sCDη-hL2GρscRqtαR2sCDDL2ηcoth⁡⁡R2sCDDL2η-1.


## 3.1. The Two-Phase Flow Governing Equations in Fracture System

Governing flow equations of two-phase flow in fractures under three-dimensional Cartesian coordinates can be written as follows.Gas flow equation in fractures is(1)∇ρgKKfrgμg∇pf-qm=∂ϕfSgρg∂t.Liquid flow equation in fractures is(2)∇ρlKKfrlμl∇pf=∂ϕfSlρl∂t.The transient diffusion rate of the shale gas from the shale matrixqm can be expressed as follows:(3)qm=G∂cm∂rm.The two-phase flow pseudopressure is employed and defined as follows:(4)mp=∫0pfρgKfrgμg+ρlKfrlμldpf.With (1), (2), and (4), the two-phase flow governing equation can be written as follows:(5)∇2m-qm=ϕfK∂Sgρg+Slρl∂t.Employing Taylor expansion to simplify (5), we can obtain the following equation:(6)∇2m-qm=η∂m∂t,where (7)η=ϕfCtKρgKfrg/μg+ρlKfrl/μl.Employing source function into (6) yields(8)∇2m+1K-qm+qgpscRTsc+qlρl∇2m+1KqgpscRTsc·δx-xi,y-yi,z-zi,t-τ=η∂m∂t.To simplify the mathematical model, dimensionless parameters are defined as follows.The dimensionless two-phase pseudopressure is(9)mD=Khqtmi-m.The dimensionless wellbore storage coefficient is(10)CD=CKhL2ηρgKfrg/μg+ρlKfrl/μl.The dimensionless time is(11)tD=CDtL2η.The dimensionless coordinate value ofx is(12)xD=xL.The dimensionless coordinate value ofy is(13)yD=yL.The dimensionless coordinate value ofz is(14)zD=zL.The dimensionless relative mass concentration is(15)cmD=1ρscci-cm.The dimensionless radial distance in matrix system is(16)rmD=rmR.Using the dimensionless parameters, (8) can be written as follows:(17)∇2mD-hL2qtGρscR∂cmD∂rmD+qgpscRTsc+qlρl∇2mD-hL2qtGρscR·δxD-xiD,yD-yiD,zD-ziD,tD-τD=CDη∂mD∂tD.The two-phase flow mathematical model in Laplace space is obtained by taking the Laplace transformation [28] of (17). The two-phase flow mathematical model is as follows:(18)∇2mD¯-hL2qtqgpscRTsc+qlρlδxD-xiD,yD-yiD,zD-ziD=sCDηmD¯+hL2GρscRqt∂cmD¯∂rmD.With the condition of sealed boundary, the dimensionless form of the initial condition is(19)lim|tD|→0mD¯=0.The dimensionless forms of the outer boundary conditions are(20)lim|xD|→0∂mD¯∂xD=0,lim|xD|→a/L∂mD¯∂xD=0,lim|yD|→0∂mD¯∂yD=0,lim|yD|→b/L∂mD¯∂yD=0,lim|zD|→0∂mD¯∂zD=0,lim|zD|→h/L∂mD¯∂zD=0.


## 3.2. The Transient Diffusion Rate Equations in Matrix System

With Fick’s second law and the assumption of spherical matrix blocks, the transient diffusion rate equation can be written as follows [29]:(21)1rm2∂∂rmDrm2∂cm∂rm=∂cm∂t.Using dimensionless parameters and Laplace transform to rewrite (21),(22)∂2cmD¯∂rmD2+2rmD∂cmD¯∂rmD=R2uCDDL2ηcmD¯.For transient diffusion case, the center of the spherical matrix block can be treated as a no-flow boundary, so the inner boundary condition can be written as follows:(23)∂cmD∂rmDrmD=0=0.The concentration on the external surface of matrix block is in equilibrium with the pressure in the fracture. The outer boundary condition can be written as follows:(24)cmDrmD,tDrmD=1=cmDm.The general solution of (22) can be expressed as follows:(25)cmD¯=ArmDeWrmD+BrmDe-WrmD.With the boundary condition and general solution expression, the accurate solution of the dimensionless mathematical model of matrix system for transient diffusion case can be obtained as follows:(26)cmD¯=cmDm¯2shR2sCD/DL2ηrmDcmD¯=·eR2sCD/DL2ηrmD-e-R2sCD/DL2ηrmD.With the definition of hyperbolic functions, (26) can be rewritten as(27)cmD¯=cmDm¯shR2sCD/DL2ηrmDshR2sCDDL2ηrmD,where cmD(m)¯ is the dimensionless equilibrium volumetric gas concentration. Employing Langmuir isotherm [30], there is(28)cmDm¯=-VLmLmi-mfmL+mfmL+mi¯cmDm¯=-qtKhVLmLmL+mfmL+mimD¯.Equation (28) can be rewritten as(29)cmDm¯=-qtKhVLmLmL+mfmL+mimD¯=-αmD¯,(30)α=qtKhVLmLmL+mfmL+mi.Inserting (29) into (27) and getting derivative with respect to rmD, the solution of the matrix model in the Laplace domain is as follows:(31)∂cmD¯∂rmD=-αmD¯R2sCDDL2ηcoth⁡R2sCDDL2η-1.


## 3.3. The Two-Phase Flow Mathematical Model for Multistage Fractured Horizontal Well in SGR

Inserting (31) into (18) results in(32)∇2mD¯-hL2qtqgpscRTsc+qlρlδ=sCDη-hL2GρscRqtαR2sCDDL2ηcoth⁡⁡R2sCDDL2η-1·mD¯.To simplify, (32) can be rewritten as(33)∇2mD¯-hL2qtqgpscRTsc+qlρlδ=fsmD¯,(34)fs=sCDη-hL2GρscRqtαR2sCDDL2ηcoth⁡⁡R2sCDDL2η-1.


## 4. Mathematical Model Solution

### 4.1. Eigenvalue Method and Orthogonal Transformation

Employing three-dimensional eigenvalue method and orthogonal transformation to solve the accurate solution of the mathematical model,(35)∂2ExD,yD,zD∂xD2+∂2ExD,yD,zD∂yD2+∂2ExD,yD,zD∂zD2=-λxyzExD,yD,zD,∂ExD,yD,zD∂xDxD=0=0,∂ExD,yD,zD∂xDxD=a/L=0,∂ExD,yD,zD∂yDyD=0=0,∂ExD,yD,zD∂yDyD=b/L=0,∂ExD,yD,zD∂zDzD=0=0,∂ExD,yD,zD∂zDzD=h/L=0,where mD¯(xD,yD,zD)=X(x)Y(y)Z(z) is the expression of dimensionless two-phase pseudopressure in three-dimensional eigenvalue space, and we can transform the three-dimensional problem to three one-dimensional problems:(36)∂2ExxD∂xD2=-λxExxD,0<xD<aL,∂ExxD∂xDxD=0=0,∂ExxD∂xDxD=a/L=0.With (36), eigenvalue in X direction can be obtained as λ=(βπL/a)2, and eigenfunction in X direction is Ex(xD)=cos⁡(βπL/axD), β=1,2,3,….In the same way,Y direction’s eigenvalue can be obtained as λ=(γπL/b)2, and eigenfunction in Y direction is EyyD=cos⁡γπL/byD, γ=1,2,3,….Z direction’s eigenvalue can be obtained as λ=(nπL/h)2, and eigenfunction in Z direction is EzzD=cos⁡nπL/hzD, n=1,2,3,….Use the one-dimensional eigenvalues and eigenfunctions to obtain three-dimensional eigenvalue and eigenfunction as follows:(37)λxyz=βπLa2+γπLb2+nπLh2(38)ExyzxD,yD,zD=cos⁡βπLaxDcos⁡γπLbyDExyzxD,yD,zD=·cos⁡nπLhzD,β,γ,n=1,2,3,….The eigenfunction in (38) is a complete orthonormal system in three-dimensional space; the completeness orthogonality equation is as follows:(39)F=∫-∞+∞∫-∞+∞∫0h/LExD,yD,zDF=∫∫∫∫∫·E′xD,yD,zDdxDdyDdzD.Take the eigenvalues and eigenfunctions into (38) to rewrite the completeness orthogonality equation as follows:(40)16Aββ′Bββ′Cγγ′Dγγ′Enn′Fnn′,where(41)Aββ′=0,β≠β′1,β=β′,Bββ′=aL,β=0a2L,β=1,2,3,…,Cγγ′=0,γ≠γ′1,γ=γ′,Dγγ′=bL,γ=0b2L,γ=1,2,3,…,Enn′=0,n≠n′1,n=n′,Fnn′=hL,n=0h2L,n=1,2,3,….EmployingF[W(x)]=∫abs(x)K(x,μ)W(x)dx to define orthogonal transformation, then (33) can be solved as follows:(42)FmD¯=hL2/qtqgpsc/RTsc+qlρlExi,yi,zifs+βπL/a2+γπL/b2+nπL/h2.With the definition of orthogonal inverse transformation and norm, there are equations as follows.Orthogonal inverse transformation is as follows:(43)Wx=∑n=1∞FWnEn2Enx.Norm is as follows:(44)En2=∫absxEn2xdx.The value of norm can be obtained as follows:(45)Enx2=12Bββ′,Eny2=12Dγγ′,Enz2=12Fnn′,En2=16Bββ′Dγγ′Fnn′.Taking (43) and (44) to transform (42), the point source expression can obtained as follows:(46)F-1mD¯=∑β=0∞∑γ=0∞∑n=0∞6hL2qgpscRTsc+qlρlExi,yi,zi+nπLh2-1∑∑∑·qtBββ′Dγγ′Fnn′βπLa2∑∑∑·fs+βπLa2+γπLb2∑∑∑qtBββ′Dγγ′6hL2qgpsc/RTsc+qlρlExi,yi,ziβπLa2+nπLh2-1ExD,yD,zD.


### 4.2. Model Solution

Based on the assumptions, all the fractures fully penetrated the SGR and are infinitely conductive. According to the Laplace transform,(47)L∫0tftdt=1sfz¯.There are transform equations as follows:(48)mDj¯=∫0tF-1mD¯dt.TheF-1[mD¯] is the point source solution of the fractures, with the property of point source function to rewrite (48), and the result can be obtained as follows:(49)mDj¯=qi¯∑β=0∞∑γ=0∞∑n=0∞6hL2Exi,yi,zi+nπLh2-1=qi¯∑∑∑·sBββ′Dγγ′Fnn′βπLa2=qi¯∑∑∑·fs+βπLa2+γπLb2=qi¯∑∑∑βπLa2+nπLh2-1Exj,yj,zj,where (50)qi¯=qgpsc/R+qwρwqt.Two-phase dimensionless characteristic production rate and other two parameters can be defined as follows:(51)qD¯=qi¯L2Δyfi,qDi=1tD∫0tDL-1[qD¯]dtD,qDi′=∂qDi∂ln⁡tD.Thus, (49) can be rewritten as follows:(52)mDj¯=qD¯∑β=0∞∑γ=0∞∑n=0∞6hExi,yi,zi+nπLh2-1=qD¯∑∑∑·sBββ′Dγγ′Fnn′βπLa2=qD¯∑∑∑·fs+βπLa2+γπLb2=qD¯∑∑∑βπLa2+nπLh2-1Exj,yj,zj.Integrating (52) with respect to z over the interval 0 to h results in(53)mDjZj¯=hqD¯Uxj,yj,zj,where(54)Uxj,yj,zj=∑β=0∞∑γ=0∞∑n=0∞6hcos⁡βπLaxjcos⁡γπLbyj+nπLh2-1=∑∑∑·Bββ′Dγγ′Fnn′βπLa2=∑∑∑·fs+βπLa2+γπLb2=∑∑∑6hcos⁡βπL/axjcos⁡γπL/byjβπLa2+nπLh2-1Exi,yi,zi.As we can see in Figure2, every fracture is divided into n equal segments in Y direction. Each discrete segment has a midpoint, (xn,yn,j+1) is the coordinate of one discrete segment, and yn,j and yn,j+1 are the start point and the end point of one discrete segment.Figure 2
Schematic of multiple transverse fractures along a horizontal well.Integrating (53) with respect to Y over the interval yn,j to yn,j+1, the result is as follows:(55)mDijxj,yj,zj¯=qD¯∫yjyj+1hUxj,yj,zjdyj.The dimensionless pseudopressure response fornf fractures can be expressed as follows:(56)mDij¯=∑i=1nf∑j=12nqD¯∫yjyj+1hUxj,yj,zjdyj.Fluid flow in wellbore and fractures is under the assumption of infinite conductivity, so the pressure at each point within the fractures and the horizontal wellbore is identical to bottom hole pressure, and there is an expression as follows:(57)mDij¯=mwD¯.With the assumption of constant flow rate, there is an expression as follows:(58)∑i=1nf∑j=12nqD¯=1s.With (58) and (56), we can establish a matrix equation group. Employing Duhamel’s theorem, the accurate solution incorporating the influence of skin factor can be obtained as follows:(59)mwDCD,SK,s¯=mwD¯+SK/2πsLD1+sCDSK/2πLD+s2CDmwD¯.The solution in Laplace space to the two-phase transient flow model at constant wellbore pressure production can be expressed as(60)qD¯=1s2mwDCD,SK,s¯.


## 4.1. Eigenvalue Method and Orthogonal Transformation

Employing three-dimensional eigenvalue method and orthogonal transformation to solve the accurate solution of the mathematical model,(35)∂2ExD,yD,zD∂xD2+∂2ExD,yD,zD∂yD2+∂2ExD,yD,zD∂zD2=-λxyzExD,yD,zD,∂ExD,yD,zD∂xDxD=0=0,∂ExD,yD,zD∂xDxD=a/L=0,∂ExD,yD,zD∂yDyD=0=0,∂ExD,yD,zD∂yDyD=b/L=0,∂ExD,yD,zD∂zDzD=0=0,∂ExD,yD,zD∂zDzD=h/L=0,where mD¯(xD,yD,zD)=X(x)Y(y)Z(z) is the expression of dimensionless two-phase pseudopressure in three-dimensional eigenvalue space, and we can transform the three-dimensional problem to three one-dimensional problems:(36)∂2ExxD∂xD2=-λxExxD,0<xD<aL,∂ExxD∂xDxD=0=0,∂ExxD∂xDxD=a/L=0.With (36), eigenvalue in X direction can be obtained as λ=(βπL/a)2, and eigenfunction in X direction is Ex(xD)=cos⁡(βπL/axD), β=1,2,3,….In the same way,Y direction’s eigenvalue can be obtained as λ=(γπL/b)2, and eigenfunction in Y direction is EyyD=cos⁡γπL/byD, γ=1,2,3,….Z direction’s eigenvalue can be obtained as λ=(nπL/h)2, and eigenfunction in Z direction is EzzD=cos⁡nπL/hzD, n=1,2,3,….Use the one-dimensional eigenvalues and eigenfunctions to obtain three-dimensional eigenvalue and eigenfunction as follows:(37)λxyz=βπLa2+γπLb2+nπLh2(38)ExyzxD,yD,zD=cos⁡βπLaxDcos⁡γπLbyDExyzxD,yD,zD=·cos⁡nπLhzD,β,γ,n=1,2,3,….The eigenfunction in (38) is a complete orthonormal system in three-dimensional space; the completeness orthogonality equation is as follows:(39)F=∫-∞+∞∫-∞+∞∫0h/LExD,yD,zDF=∫∫∫∫∫·E′xD,yD,zDdxDdyDdzD.Take the eigenvalues and eigenfunctions into (38) to rewrite the completeness orthogonality equation as follows:(40)16Aββ′Bββ′Cγγ′Dγγ′Enn′Fnn′,where(41)Aββ′=0,β≠β′1,β=β′,Bββ′=aL,β=0a2L,β=1,2,3,…,Cγγ′=0,γ≠γ′1,γ=γ′,Dγγ′=bL,γ=0b2L,γ=1,2,3,…,Enn′=0,n≠n′1,n=n′,Fnn′=hL,n=0h2L,n=1,2,3,….EmployingF[W(x)]=∫abs(x)K(x,μ)W(x)dx to define orthogonal transformation, then (33) can be solved as follows:(42)FmD¯=hL2/qtqgpsc/RTsc+qlρlExi,yi,zifs+βπL/a2+γπL/b2+nπL/h2.With the definition of orthogonal inverse transformation and norm, there are equations as follows.Orthogonal inverse transformation is as follows:(43)Wx=∑n=1∞FWnEn2Enx.Norm is as follows:(44)En2=∫absxEn2xdx.The value of norm can be obtained as follows:(45)Enx2=12Bββ′,Eny2=12Dγγ′,Enz2=12Fnn′,En2=16Bββ′Dγγ′Fnn′.Taking (43) and (44) to transform (42), the point source expression can obtained as follows:(46)F-1mD¯=∑β=0∞∑γ=0∞∑n=0∞6hL2qgpscRTsc+qlρlExi,yi,zi+nπLh2-1∑∑∑·qtBββ′Dγγ′Fnn′βπLa2∑∑∑·fs+βπLa2+γπLb2∑∑∑qtBββ′Dγγ′6hL2qgpsc/RTsc+qlρlExi,yi,ziβπLa2+nπLh2-1ExD,yD,zD.


## 4.2. Model Solution

Based on the assumptions, all the fractures fully penetrated the SGR and are infinitely conductive. According to the Laplace transform,(47)L∫0tftdt=1sfz¯.There are transform equations as follows:(48)mDj¯=∫0tF-1mD¯dt.TheF-1[mD¯] is the point source solution of the fractures, with the property of point source function to rewrite (48), and the result can be obtained as follows:(49)mDj¯=qi¯∑β=0∞∑γ=0∞∑n=0∞6hL2Exi,yi,zi+nπLh2-1=qi¯∑∑∑·sBββ′Dγγ′Fnn′βπLa2=qi¯∑∑∑·fs+βπLa2+γπLb2=qi¯∑∑∑βπLa2+nπLh2-1Exj,yj,zj,where (50)qi¯=qgpsc/R+qwρwqt.Two-phase dimensionless characteristic production rate and other two parameters can be defined as follows:(51)qD¯=qi¯L2Δyfi,qDi=1tD∫0tDL-1[qD¯]dtD,qDi′=∂qDi∂ln⁡tD.Thus, (49) can be rewritten as follows:(52)mDj¯=qD¯∑β=0∞∑γ=0∞∑n=0∞6hExi,yi,zi+nπLh2-1=qD¯∑∑∑·sBββ′Dγγ′Fnn′βπLa2=qD¯∑∑∑·fs+βπLa2+γπLb2=qD¯∑∑∑βπLa2+nπLh2-1Exj,yj,zj.Integrating (52) with respect to z over the interval 0 to h results in(53)mDjZj¯=hqD¯Uxj,yj,zj,where(54)Uxj,yj,zj=∑β=0∞∑γ=0∞∑n=0∞6hcos⁡βπLaxjcos⁡γπLbyj+nπLh2-1=∑∑∑·Bββ′Dγγ′Fnn′βπLa2=∑∑∑·fs+βπLa2+γπLb2=∑∑∑6hcos⁡βπL/axjcos⁡γπL/byjβπLa2+nπLh2-1Exi,yi,zi.As we can see in Figure2, every fracture is divided into n equal segments in Y direction. Each discrete segment has a midpoint, (xn,yn,j+1) is the coordinate of one discrete segment, and yn,j and yn,j+1 are the start point and the end point of one discrete segment.Figure 2
Schematic of multiple transverse fractures along a horizontal well.Integrating (53) with respect to Y over the interval yn,j to yn,j+1, the result is as follows:(55)mDijxj,yj,zj¯=qD¯∫yjyj+1hUxj,yj,zjdyj.The dimensionless pseudopressure response fornf fractures can be expressed as follows:(56)mDij¯=∑i=1nf∑j=12nqD¯∫yjyj+1hUxj,yj,zjdyj.Fluid flow in wellbore and fractures is under the assumption of infinite conductivity, so the pressure at each point within the fractures and the horizontal wellbore is identical to bottom hole pressure, and there is an expression as follows:(57)mDij¯=mwD¯.With the assumption of constant flow rate, there is an expression as follows:(58)∑i=1nf∑j=12nqD¯=1s.With (58) and (56), we can establish a matrix equation group. Employing Duhamel’s theorem, the accurate solution incorporating the influence of skin factor can be obtained as follows:(59)mwDCD,SK,s¯=mwD¯+SK/2πsLD1+sCDSK/2πLD+s2CDmwD¯.The solution in Laplace space to the two-phase transient flow model at constant wellbore pressure production can be expressed as(60)qD¯=1s2mwDCD,SK,s¯.


## 5. Analysis of Type Curves Characteristics

### 5.1. The Type Curves of Pseudopressure Response and Production Decline

The type curves of a two-phase multifractured horizontal well in SGR can be plotted by Stehfest [31] numerical invention algorithm. Intuitively, the type curves graphically show the process and characteristics of fluid flow in reservoirs.In the condition of sealed boundary, the two-phase pseudopressure-transient pressure response of multistage fractured horizontal well in SGR, which has six flow regimes, can be clearly shown in Figure3. Regime 1 is the pure wellbore storage regime; the pseudopressure curve and its derivative curve have the same trend as conventional gas reservoir. Regime 2 is hump transition regime, which is caused by fluids difference in two-phase flow; the shape of the derivative curve looks like a concave. Regime 3 is early linear flow regime; the flow energy in this regime is provided by free gas flow. Formation pressure is still over desorption critical pressure, so the impact of desorbed gas cannot be seen in this period. Regime 4 is early radial flow regime, which behaves as a horizontal line in type curve. The clarity of this period primarily depends on the fracture length and spacing. When facture spacing is small enough and fracture length is long enough, this period may not be so obvious in type curve. Regime 5 is desorbed gas transient diffusion flow regime, and the derivative curve looks like an “arch,” which is caused by desorption. It is obviously influenced by absorption index. Regime 6 is late-time compound flow regime, and because of sealed boundary condition, it behaves like an up-warp in the type curve.Figure 3
Type curve of pressure transient response for two-phase flow in multistage hydraulic fractured horizontal well in SGR.As we can see in Figure4, the dimensionless production decline type curve is plotted by using the Stehfest numerical invention algorithm and Blasingame production decline analysis method. The dimensionless production decline type curve can be divided into six flow stages:Stage 1: wellbore storage which has the same characteristics as the conventional reservoir;Stage 2: formation linear flow period, which is characterized by a horizontal line on both dimensionless production integral and its derivative curves;Stage 3: interporosity flow period, which is characterized by a V shape in the derivative curve, and in this period free gas has exhausted, and desorbed gas from the shale matrix surface diffuses to factures by pressure difference;Stage 4: fracture linear flow period, which is characterized by a drop line in type curves;Stage 5: compound radial flow period, which is characterized by a drop line in dimensionless production integral curve and a rising line in dimensionless production derivative curve;Stage 6: boundary control flow period, which is characterized by a rapid drop line on both dimensionless production integral and its derivative curves for the reason of sealed boundary condition.Figure 4
Type curve of production decline for two-phase flow in multistage hydraulic fractured horizontal well in SGR.


### 5.2. Sensitivity Analysis of Production Decline

Figures5–9 show the type curves of production decline analysis for a two-phase multifractured horizontal well in SGR. These figures are the typical response of two-phase multifractured horizontal well percolation in SGR at constant wellbore pressure production.Figure 5
Effect of initial saturation on type curves.Figure 6
Effect of skin factor on type curves.Figure 7
Effect of absorption index on type curves.Figure 8
Effect of different fracture stages on type curves.Figure 9
Effect of horizontal well lateral length on type curves.Figure5 shows the production decline type curve characteristics affected by the initial saturation. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the increase of gas initial saturation.Figure6 shows the production decline type curve characteristics affected by the skin factor. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the increase of skin factor. This characteristic can be obviously seen in the early two stages: wellbore storage period and formation linear flow period.Figure7 shows the production decline type curve characteristics affected by the absorption index. It can be seen from the figure that the interporosity flow period may appear earlier with the increase of skin factor, and the fracture linear flow period will last even longer. The effort of absorption index cannot be observed in early product period and later product period.Figure8 shows the production decline type curve characteristics affected by the different fracture stages. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the decrease of fracture stages. More fractures would make flow condition near the horizontal well better, but if there are relative large fracture stages, the increase rate of the type curves may not be so obvious.Figure9 shows the production decline type curve characteristics affected by the horizontal well lateral length. It can be seen from the figure that the effort of horizontal well lateral length can be obviously seen in the early two stages, especially when the horizontal well lateral length is relatively long. If the horizontal well lateral length is long enough, the dimensionless production rate derivative curve may go down later and be smooth in the early flow period.


## 5.1. The Type Curves of Pseudopressure Response and Production Decline

The type curves of a two-phase multifractured horizontal well in SGR can be plotted by Stehfest [31] numerical invention algorithm. Intuitively, the type curves graphically show the process and characteristics of fluid flow in reservoirs.In the condition of sealed boundary, the two-phase pseudopressure-transient pressure response of multistage fractured horizontal well in SGR, which has six flow regimes, can be clearly shown in Figure3. Regime 1 is the pure wellbore storage regime; the pseudopressure curve and its derivative curve have the same trend as conventional gas reservoir. Regime 2 is hump transition regime, which is caused by fluids difference in two-phase flow; the shape of the derivative curve looks like a concave. Regime 3 is early linear flow regime; the flow energy in this regime is provided by free gas flow. Formation pressure is still over desorption critical pressure, so the impact of desorbed gas cannot be seen in this period. Regime 4 is early radial flow regime, which behaves as a horizontal line in type curve. The clarity of this period primarily depends on the fracture length and spacing. When facture spacing is small enough and fracture length is long enough, this period may not be so obvious in type curve. Regime 5 is desorbed gas transient diffusion flow regime, and the derivative curve looks like an “arch,” which is caused by desorption. It is obviously influenced by absorption index. Regime 6 is late-time compound flow regime, and because of sealed boundary condition, it behaves like an up-warp in the type curve.Figure 3
Type curve of pressure transient response for two-phase flow in multistage hydraulic fractured horizontal well in SGR.As we can see in Figure4, the dimensionless production decline type curve is plotted by using the Stehfest numerical invention algorithm and Blasingame production decline analysis method. The dimensionless production decline type curve can be divided into six flow stages:Stage 1: wellbore storage which has the same characteristics as the conventional reservoir;Stage 2: formation linear flow period, which is characterized by a horizontal line on both dimensionless production integral and its derivative curves;Stage 3: interporosity flow period, which is characterized by a V shape in the derivative curve, and in this period free gas has exhausted, and desorbed gas from the shale matrix surface diffuses to factures by pressure difference;Stage 4: fracture linear flow period, which is characterized by a drop line in type curves;Stage 5: compound radial flow period, which is characterized by a drop line in dimensionless production integral curve and a rising line in dimensionless production derivative curve;Stage 6: boundary control flow period, which is characterized by a rapid drop line on both dimensionless production integral and its derivative curves for the reason of sealed boundary condition.Figure 4
Type curve of production decline for two-phase flow in multistage hydraulic fractured horizontal well in SGR.


## 5.2. Sensitivity Analysis of Production Decline

Figures5–9 show the type curves of production decline analysis for a two-phase multifractured horizontal well in SGR. These figures are the typical response of two-phase multifractured horizontal well percolation in SGR at constant wellbore pressure production.Figure 5
Effect of initial saturation on type curves.Figure 6
Effect of skin factor on type curves.Figure 7
Effect of absorption index on type curves.Figure 8
Effect of different fracture stages on type curves.Figure 9
Effect of horizontal well lateral length on type curves.Figure5 shows the production decline type curve characteristics affected by the initial saturation. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the increase of gas initial saturation.Figure6 shows the production decline type curve characteristics affected by the skin factor. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the increase of skin factor. This characteristic can be obviously seen in the early two stages: wellbore storage period and formation linear flow period.Figure7 shows the production decline type curve characteristics affected by the absorption index. It can be seen from the figure that the interporosity flow period may appear earlier with the increase of skin factor, and the fracture linear flow period will last even longer. The effort of absorption index cannot be observed in early product period and later product period.Figure8 shows the production decline type curve characteristics affected by the different fracture stages. It can be seen from the figure that the curves of the dimensionless production rate integral and its derivative curves go down with the decrease of fracture stages. More fractures would make flow condition near the horizontal well better, but if there are relative large fracture stages, the increase rate of the type curves may not be so obvious.Figure9 shows the production decline type curve characteristics affected by the horizontal well lateral length. It can be seen from the figure that the effort of horizontal well lateral length can be obviously seen in the early two stages, especially when the horizontal well lateral length is relatively long. If the horizontal well lateral length is long enough, the dimensionless production rate derivative curve may go down later and be smooth in the early flow period.


## 6. Conclusions

The transient flow model of two-phase multifractured horizontal well and desorbed gas transient diffusion in SGR is established and solved, the type curves of pressure response and production decline are illustrated, and the production decline performance under the effects of hydraulic fractures and SGR properties are discussed. The following conclusions were obtained.(1)
The two-phase multifractured horizontal well model can be solved using the eigenvalue method and orthogonal transformation. The desorbed gas transient diffusion rate can be solved using Langmuir’s equilibrium sorption isotherm and Fick’s second law.(2)
The type curves of production decline are dominated by initial saturation, skin factor, and fracture stages. The curves of the dimensionless production rate integral and its derivative curves go down with the increase of gas initial saturation. The skin factor has primarily effect on the early two stages; the larger the skin factor is, the lower the derivative curves are. More fracture stages can improve the flow condition near a horizontal well, and the number of stages indicate the value of fracture spacing; larger fracture spacing increases the duration of fluid flow around individual fractures, and the more the fracture stages are, the more the dimensionless production rate is.(3)
The absorption index has a primary effect on the starting time of the interporosity flow period and the depth of the trough. The trough in the type curves is an indication of desorbed gas diffusion in SGR; the more the absorption index is, the later the decline appears, and, in addition, the trough becomes more significant with larger value of absorption index.(4)
The horizontal well lateral length has a primary effect on the early two flow stages but does not affect the intermediate and later flow behavior; the smaller the horizontal well lateral length is, the faster the decline appears.



---

*Source: 212103-2015-09-27.xml*