# Pressure Transient Analysis and Flux Distribution for Multistage Fractured Horizontal Wells in Triple-Porosity Reservoir Media with Consideration of Stress-Sensitivity Effect

Authors: Jingjing Guo; Haitao Wang; Liehui Zhang; Chengyong Li
Journal: Journal of Chemistry
 (2015)

---

## Abstract

Triple-porosity model is usually adopted to describe reservoirs with multiscaled pore spaces, including matrix pores, natural fractures, and vugs. Multiple fractures created by hydraulic fracturing can effectively improve the connectivity between existing natural fractures and thus increase well deliverability. However, little work has been done on pressure transient behavior of multistage fractured horizontal wells in triple-porosity reservoirs. Based on source/sink function method, this paper presents a triple-porosity model to investigate the transient pressure dynamics and flux distribution for multistage fractured horizontal wells in fractured-vuggy reservoirs with consideration of stress-dependent natural fracture permeability. The model is semianalytically solved by discretizing hydraulic fractures and Pedrosa’s transformation, perturbation theory, and integration transformation method. Type curves of transient pressure dynamics are generated, and flux distribution among hydraulic fractures for a fractured horizontal well with constant production rate is also discussed. Parametric study shows that major influential parameters on transient pressure responses are parameters pertinent to reservoir properties, interporosity mass transfer, and hydraulic fractures. Analysis of flux distribution indicates that flux density gradually increases from the horizontal wellbore to fracture tips, and the flux contribution of outermost fractures is higher than that of inner fractures. The model can also be extended to optimize hydraulic fracture parameters.


---

## Body

## 1. Introduction

The concept of dual-porosity media was first proposed by Barenblatt et al. [1] to describe naturally fractured media. Warren and Root [2] then extended this model to analyze pressure transient dynamics for vertical wells in naturally fractured oil reservoirs. Recently, Cai et al. [3, 4] analyzed the imbibition mechanism in fractured-porous media based on fractal geometry. Although the dual-porosity model can describe most naturally fractured reservoirs, there are still some drawbacks when applying the dual-porosity model to describe reservoirs with multiple pore media, such as reservoirs with natural fractures and vugs (different from matrix pores).In order to better describe reservoir heterogeneity in this type of naturally fractured-vuggy reservoirs, the concept of triple-porosity model was proposed. The first triple-porosity model in the literature was proposed by Abadassah and Ershaghi [5], in which the reservoir was represented by a combination of two matrix systems with different properties and uniformly distributed natural fractures. The model was then further extended by Al-Ghamdi and Ershaghi [6] to represent reservoirs with two fracture systems (microfractures and macrofractures) and one matrix system. After that, over the past several decades, numerous efforts have been made to investigate transient pressure behavior in so-called triple-porosity reservoirs. However, almost all the research focuses on vertical wells [7–19] or horizontal wells [20–22] in triple-porosity reservoirs, and few studies have been done on complex well-reservoir configuration, such as multistage fractured horizontal well in triple-porosity reservoirs.Field practices have proven that horizontal wellbore combining multistage hydraulic fracturing can effectively enhance the permeability in the vicinity of horizontal wellbore. Almost all the triple-porosity models [23, 24] for fractured horizontal wells assume linear flow, which cannot completely characterize the pressure responses during all flowing periods (such as pseudoradial flow period). On the other hand, in most existing models for fractured horizontal wells in triple-porosity reservoirs three porous media denote hydraulic fractures, natural fractures, and matrix pores, which means the reservoir itself is actually considered as a dual-porosity medium.This study develops a semianalytical model based on source/sink function theory for analyzing transient pressure responses and flux distribution of naturally fractured-vuggy reservoirs, in which the reservoir itself is conceptualized as triple-porosity media and the interaction between hydraulic fractures and reservoir is considered. The model presented here can completely reflect transient pressure characteristics during all the possible flowing period as well as flux distribution among multiple fractures. In addition, this model also takes into account the stress-sensitivity of permeability caused by closure of natural fractures during production.


## 2. Physical Model of a Multistage Fractured Horizontal Well in a Triple-Porosity Reservoir

As shown in Figure1, a horizontal well with multiple hydraulic fractures, producing at a constant rate (q) in a reservoir with multiscale storage spaces, is considered here. Both the cap rock and underlying rock formation of the reservoir are usually shales with extremely low permeability, which can be conceived as no-flow boundaries. The reservoir thickness and oil viscosity are assumed to be h and μ, separately. Other basic assumptions involved in the derivation of the triple-porosity model are as follows:(1)
The reservoir is assumed to be composed of three contiguous porous media, which are matrix, natural fractures, and vugs (representing larger voids in the reservoir which are not part of natural fractures).(2)
Natural fractures are assumed to be directly connected to hydraulic fractures or the horizontal wellbore, while the fluid stored in matrix pores or vugs does not have direct access to hydraulic fractures or the horizontal wellbore. Figure2 presents a simple illustration of the fluid flowing path in triple-porosity reservoirs.(3)
The permeability of natural fractures is considered stress-sensitive to incorporate the effect of partial or complete closure of natural fractures during production.(4)
Both oil and rock are considered slightly compressible, and the reservoir has uniform pressure distribution at time zero.(5)
Consider single-phase isothermal flow and negligible gravity and capillary effects.Figure 1
Schematic of a multistage fractured horizontal well in a triple-porosity reservoir.Figure 2
Illustration of fluid flowing paths in a triple-porosity reservoir.To obtain the pressure responses caused by the well-reservoir configuration shown in Figure1, a mathematical model together with corresponding boundary and initial conditions is required. In the problem addressed here, the multiple fractures inner boundary condition is so complex that it is rather difficult to directly describe it with a mathematical equation. In order to solve this problem, source function theory, a powerful tool to solve complex transient flow problems in reservoirs, is adopted in this work. We first studied the transient pressure responses caused by a continuous line-sink in triple-porosity reservoirs and then adopted the line-sink solution with superposition principle to obtain the pressure transient dynamics and flux distribution of multistage fractured horizontal wells in triple-porosity reservoirs.


## 3. Line-Sink Model in Stress-Sensitive Triple-Porosity Reservoirs

### 3.1. Physical Model for a Line-Sink in Triple-Porosity Reservoirs

As shown in Figure3, a continuous line-sink in a triple-porosity reservoir with stress-dependent fracture permeability is considered in this section. The strength of the line-sink is represented by q^(t).Figure 3
A line-sink in a stress-sensitive triple-porosity reservoir.


### 3.2. Mathematical Model for a Line-Sink in Triple-Porosity Reservoirs

(1) Governing Equations. Imposing mass-balance equation over a control volume in the triple-porosity reservoir and substituting equation of state and equation of motion, we can get the governing equations for matrix, natural fracture, and vug systems.For natural fracture system with stress-dependent permeability, the governing equation is as follows (i.e., (A.26b) in Appendix A):(1)∂2pf∂r2+1r∂pf∂r+γ∂pf∂r2=e-γpf-piϕf0μCftkfi∂pf∂t+ϕm0μCmtkfi∂pm∂t+ϕv0μCvtkfi∂pv∂t.For matrix system, the governing equation is as follows (i.e., (A.30) in Appendix A):(2)ϕm0Cmt∂pm∂t+α1kmμpm-pf=0.For vug system, the governing equation is as follows (i.e., (A.34) in Appendix A):(3)ϕv0Cvt∂pv∂t+α2kvμpv-pf=0.(2) Boundary Conditions. According to the derivation in Appendix B, the inner boundary condition of the line-sink model is (i.e., (B.5) in Appendix B)(4)limε→0⁡e-γpi-pfr∂pf∂rr=ε=q^tμB2πkfih.Assuming infinitely large reservoirs, the outer boundary condition can be expressed as(5)pfr,tr→∞=pi.(3) Initial Condition. The initial pressure distribution in the reservoir is assumed to be uniform, which is(6)pfr,tt=0=pmr,tt=0=pvr,tt=0=pi.


### 3.3. Dimensionless Form of the Line-Sink Model

For the purpose of simplifying derivation and comparing between different reservoirs, the mathematical model is derived and solved in dimensionless form in this work. The pressure drops involved in the mathematical model areΔpf=pi-pf, Δpm=pi-pm, and Δpv=pi-pv, and relevant dimensionless variables used in the mathematical model are listed in Table 1.Table 1
Definitions of dimensionless variables.
Dimensionless pressure
p
f
D
=
2
π
k
f
i
h
q
B
μ
Δ
p
f, pmD=2πkfihqBμΔpm, pvD=2πkfihqBμΔpv
Dimensionless radial distance
r
D
=
r
L, εD=εL
Dimensionless time
t
D
=
k
f
i
t
μ
ϕ
f
0
C
f
t
+
ϕ
m
0
C
m
t
+
ϕ
v
0
C
v
t
L
2
Dimensionless storativity ratio of natural fracture system
ω
f
=
ϕ
f
0
C
f
t
ϕ
f
0
C
f
t
+
ϕ
m
0
C
m
t
+
ϕ
v
0
C
v
t
Dimensionless storativity ratio of matrix system
ω
m
=
ϕ
m
0
C
m
t
ϕ
f
0
C
f
t
+
ϕ
m
0
C
m
t
+
ϕ
v
0
C
v
t
Dimensionless storativity ratio of vug system
ω
v
=
ϕ
v
0
C
v
t
ϕ
f
0
C
f
t
+
ϕ
m
0
C
m
t
+
ϕ
v
0
C
v
t
Dimensionless interporosity flow coefficient between matrix and natural fractures
λ
m
f
=
α
m
f
k
m
L
2
k
f
i
Dimensionless interporosity flow coefficient between vugs and natural fractures
λ
v
f
=
α
v
f
k
v
L
2
k
f
i
Dimensionless permeability modulus
γ
D
=
q
B
μ
2
π
k
f
i
h
γ
Dimensionless production rate of the line-sink
q
^
D
=
q
^
D
(
t
D
)
=
q
^
t
q
q
^
t
q
X
∑
X
X
X
X
Note:L is an arbitrary reference length which can be rw, Lh, and so forth.Obviously,ωf, ωm, and ωv in Table 1 satisfy the following equation:(7)ωf+ωm+ωv=1.With the definitions in Table1, the dimensionless forms of (1)~(6) can be obtained as follows:(8)∂2pfD∂rD2+1rD∂pfD∂rD-γD∂pfD∂rD2=eγDpfDωf∂pfD∂tD+ωm∂pmD∂tD+ωv∂pvD∂tD,(9)ωm∂pmD∂tD+λmfpmD-pfD=0,(10)1-ωf-ωm∂pvD∂tD+λvfpvD-pfD=0,(11)limεD→0⁡e-pfDγDrD∂pfD∂rDrD=εD=-q^D,(12)pfDrD→∞=0,(13)pfDtD=0=pmDtD=0=pvDtD=0=0.


### 3.4. Solution of the Dimensionless Line-Sink Model

#### 3.4.1. Pedrosa’s Linearization

The stress-sensitivity of natural fracture permeability makes the above seepage model strongly nonlinear which cannot be solved analytically. Here, following Pedrosa Jr. [25], we introduced the following equation to alleviate the nonlinearity of (8) and (11):(14)pfDrD,tD=-1γDln⁡1-γDUDrD,tD,where UD is an intermediate dimensionless parameter which is also a function of radial distance and time.With (14), (8)~(13) become(15)∂2UD∂rD2+1rD∂UD∂rD=11-γDUDωf∂UD∂tD+ωm∂pmD∂tD+ωv∂pvD∂tD,(16)ωm∂pmD∂tD+λmfpmD+1γDln⁡1-γDUD=0,(17)1-ωf-ωm∂pvD∂tD+λvfpvD+1γDln⁡1-γDUD=0,(18)limεD→0⁡rD∂UD∂rDrD=εD=-q^D,(19)UDrD→∞=0,(20)UDtD=0=pmDtD=0=pvDtD=0=0.


#### 3.4.2. Perturbation Technique

According to the perturbation theory, the termsUD, -(1/γD)ln⁡[1-γDUD(rD,tD)], and 1/(1-γDUD(rD,tD)) in (15) through (20) can be expanded as power series in dimensionless permeability modulus, and we can get the following equations:(21)UD=UD0+γDUD1+γD2UD2+⋯,-1γDln⁡1-γDUDrD,tD=UDrD,tD+12γDUD2rD,tD+⋯,11-γDUDrD,tD=1+γDUDrD,tD+γD2UDrD,tD+⋯Given that the dimensionless permeability modulusγD is usually small (much smaller than 1), the zero-order approximate solution can satisfy the requirements of engineering precision, so (15)~(20) become(22)∂2UD0∂rD2+1rD∂UD0∂rD=ωf∂UD0∂tD+ωm∂pmD∂tD+ωv∂pvD∂tD,ωm∂pmD∂tD+λmfpmD-UD0=0,1-ωf-ωm∂pvD∂tD+λvfpvD-UD0=0,limεD→0⁡rD∂UD0∂rDrD=εD=-q^D,UD0rD→∞=0,UD0tD=0=pmDtD=0=pvDtD=0=0.Equation (22) is the linearized line-sink model in triple-porosity reservoirs.


#### 3.4.3. Laplace Transformation

Take the following Laplace transformation:(23)UD0¯=∫0+∞UD0e-utDdtD,where u is the Laplace transformation variable.Then (22) becomes(24)d2UD0¯drD2+1rDdUD0¯drD=ωfuUD0¯+ωmupmD¯+ωvupvD¯,(25)ωmupmD¯+λmfpmD¯-UD0¯=0,(26)ωvupvD¯+λvfpvD¯-UD0¯=0,(27)limεD→0⁡rDdUD0¯drDrD=εD=-q^D¯,(28)UD0¯rD→∞=0.


#### 3.4.4. Line-Sink Solution

Equations (25) and (26) can be rewritten as(29)pmD¯=λmfωmu+λmfUD0¯,pvD¯=λvfωvu+λvfUD0¯.Substituting (29) into (24) yields(30)d2UD0¯drD2+1rDdUD0¯drD=ωf+λmfωmu+λmfωm+λvfωvu+λvfωvuUD0¯.If we define(31)gu=ωf+λmfωmu+λmfωm+λvfωvu+λvfωvu,then (30) can be expressed as follows:(32)d2UD0¯drD2+1rDdUD0¯drD=guUD0¯.The general solution of (32) can be easily obtained as follows:(33)UD0¯=c1K0gurD+c2I0gurD,where K0 and I0 are Bessel functions and c1 and c2 are unknown coefficients which can be determined by corresponding boundary conditions.If the line-sink is located at the origin of coordinates (center of the triple-porosity reservoir), then the dimensionless radial distancerD in (31) can be calculated by the following equation:(34a)rD=xD2+yD2.If the line-sink is located at(xw,yw) instead, then rD is calculated by the following equation:(34b)rD=xD-xwD2+yD-ywD2,where(35)xD=xL,(36)yD=yL,(37)xwD=xwL,(38)ywD=ywL.Substitution of (33) into the inner boundary condition (27) yields(39)limεD→0⁡gurD·-c1K1gurD+c2I1gurDrD=εD=-q^D¯.According to the properties of Bessel’s function, which arelimx→0⁡xK1(x)→1 and limx→0⁡xI1(x)→0, (39) becomes (40)c1=q^D¯.Similarly, with the outer boundary condition given in (28) and the properties of Bessel’s function, limx→∞⁡K0(x)→0 and limx→∞⁡I0(x)→∞, the coefficient in (33), c2, should satisfy the following equation:(41)c2=0.With (40) and (41), we can get the final form of (33) as follows:(42)UD0¯=q^D¯K0gurD.Equation (42) is the basic zero-order perturbation solution for a line-sink in a stress-sensitive triple-porosity reservoir. As mentioned above, zero-order perturbation solution is enough to approximate the exact solution of (15) to (20); that is, UD¯≈UD0¯. With the definition of q^D in Table 1, (42) can be expressed as(43)UD¯=q^¯qK0gurD.Equation (43) can be rewritten as(44)UD¯=q^¯qΩxD,yD;xwD,ywD,where(45)ΩxD,yD;xwD,ywD=K0gurD.Equation (44) is the basic line-sink solution for stress-sensitive triple-porosity reservoirs. With the basic line-sink solution and the superposition principle, we can obtain the pressure responses caused by multistage fractured horizontal wells in stress-sensitive triple-porosity reservoirs.


## 3.1. Physical Model for a Line-Sink in Triple-Porosity Reservoirs

As shown in Figure3, a continuous line-sink in a triple-porosity reservoir with stress-dependent fracture permeability is considered in this section. The strength of the line-sink is represented by q^(t).Figure 3
A line-sink in a stress-sensitive triple-porosity reservoir.


## 3.2. Mathematical Model for a Line-Sink in Triple-Porosity Reservoirs

(1) Governing Equations. Imposing mass-balance equation over a control volume in the triple-porosity reservoir and substituting equation of state and equation of motion, we can get the governing equations for matrix, natural fracture, and vug systems.For natural fracture system with stress-dependent permeability, the governing equation is as follows (i.e., (A.26b) in Appendix A):(1)∂2pf∂r2+1r∂pf∂r+γ∂pf∂r2=e-γpf-piϕf0μCftkfi∂pf∂t+ϕm0μCmtkfi∂pm∂t+ϕv0μCvtkfi∂pv∂t.For matrix system, the governing equation is as follows (i.e., (A.30) in Appendix A):(2)ϕm0Cmt∂pm∂t+α1kmμpm-pf=0.For vug system, the governing equation is as follows (i.e., (A.34) in Appendix A):(3)ϕv0Cvt∂pv∂t+α2kvμpv-pf=0.(2) Boundary Conditions. According to the derivation in Appendix B, the inner boundary condition of the line-sink model is (i.e., (B.5) in Appendix B)(4)limε→0⁡e-γpi-pfr∂pf∂rr=ε=q^tμB2πkfih.Assuming infinitely large reservoirs, the outer boundary condition can be expressed as(5)pfr,tr→∞=pi.(3) Initial Condition. The initial pressure distribution in the reservoir is assumed to be uniform, which is(6)pfr,tt=0=pmr,tt=0=pvr,tt=0=pi.


## 3.3. Dimensionless Form of the Line-Sink Model

For the purpose of simplifying derivation and comparing between different reservoirs, the mathematical model is derived and solved in dimensionless form in this work. The pressure drops involved in the mathematical model areΔpf=pi-pf, Δpm=pi-pm, and Δpv=pi-pv, and relevant dimensionless variables used in the mathematical model are listed in Table 1.Table 1
Definitions of dimensionless variables.
Dimensionless pressure
p
f
D
=
2
π
k
f
i
h
q
B
μ
Δ
p
f, pmD=2πkfihqBμΔpm, pvD=2πkfihqBμΔpv
Dimensionless radial distance
r
D
=
r
L, εD=εL
Dimensionless time
t
D
=
k
f
i
t
μ
ϕ
f
0
C
f
t
+
ϕ
m
0
C
m
t
+
ϕ
v
0
C
v
t
L
2
Dimensionless storativity ratio of natural fracture system
ω
f
=
ϕ
f
0
C
f
t
ϕ
f
0
C
f
t
+
ϕ
m
0
C
m
t
+
ϕ
v
0
C
v
t
Dimensionless storativity ratio of matrix system
ω
m
=
ϕ
m
0
C
m
t
ϕ
f
0
C
f
t
+
ϕ
m
0
C
m
t
+
ϕ
v
0
C
v
t
Dimensionless storativity ratio of vug system
ω
v
=
ϕ
v
0
C
v
t
ϕ
f
0
C
f
t
+
ϕ
m
0
C
m
t
+
ϕ
v
0
C
v
t
Dimensionless interporosity flow coefficient between matrix and natural fractures
λ
m
f
=
α
m
f
k
m
L
2
k
f
i
Dimensionless interporosity flow coefficient between vugs and natural fractures
λ
v
f
=
α
v
f
k
v
L
2
k
f
i
Dimensionless permeability modulus
γ
D
=
q
B
μ
2
π
k
f
i
h
γ
Dimensionless production rate of the line-sink
q
^
D
=
q
^
D
(
t
D
)
=
q
^
t
q
q
^
t
q
X
∑
X
X
X
X
Note:L is an arbitrary reference length which can be rw, Lh, and so forth.Obviously,ωf, ωm, and ωv in Table 1 satisfy the following equation:(7)ωf+ωm+ωv=1.With the definitions in Table1, the dimensionless forms of (1)~(6) can be obtained as follows:(8)∂2pfD∂rD2+1rD∂pfD∂rD-γD∂pfD∂rD2=eγDpfDωf∂pfD∂tD+ωm∂pmD∂tD+ωv∂pvD∂tD,(9)ωm∂pmD∂tD+λmfpmD-pfD=0,(10)1-ωf-ωm∂pvD∂tD+λvfpvD-pfD=0,(11)limεD→0⁡e-pfDγDrD∂pfD∂rDrD=εD=-q^D,(12)pfDrD→∞=0,(13)pfDtD=0=pmDtD=0=pvDtD=0=0.


## 3.4. Solution of the Dimensionless Line-Sink Model

### 3.4.1. Pedrosa’s Linearization

The stress-sensitivity of natural fracture permeability makes the above seepage model strongly nonlinear which cannot be solved analytically. Here, following Pedrosa Jr. [25], we introduced the following equation to alleviate the nonlinearity of (8) and (11):(14)pfDrD,tD=-1γDln⁡1-γDUDrD,tD,where UD is an intermediate dimensionless parameter which is also a function of radial distance and time.With (14), (8)~(13) become(15)∂2UD∂rD2+1rD∂UD∂rD=11-γDUDωf∂UD∂tD+ωm∂pmD∂tD+ωv∂pvD∂tD,(16)ωm∂pmD∂tD+λmfpmD+1γDln⁡1-γDUD=0,(17)1-ωf-ωm∂pvD∂tD+λvfpvD+1γDln⁡1-γDUD=0,(18)limεD→0⁡rD∂UD∂rDrD=εD=-q^D,(19)UDrD→∞=0,(20)UDtD=0=pmDtD=0=pvDtD=0=0.


### 3.4.2. Perturbation Technique

According to the perturbation theory, the termsUD, -(1/γD)ln⁡[1-γDUD(rD,tD)], and 1/(1-γDUD(rD,tD)) in (15) through (20) can be expanded as power series in dimensionless permeability modulus, and we can get the following equations:(21)UD=UD0+γDUD1+γD2UD2+⋯,-1γDln⁡1-γDUDrD,tD=UDrD,tD+12γDUD2rD,tD+⋯,11-γDUDrD,tD=1+γDUDrD,tD+γD2UDrD,tD+⋯Given that the dimensionless permeability modulusγD is usually small (much smaller than 1), the zero-order approximate solution can satisfy the requirements of engineering precision, so (15)~(20) become(22)∂2UD0∂rD2+1rD∂UD0∂rD=ωf∂UD0∂tD+ωm∂pmD∂tD+ωv∂pvD∂tD,ωm∂pmD∂tD+λmfpmD-UD0=0,1-ωf-ωm∂pvD∂tD+λvfpvD-UD0=0,limεD→0⁡rD∂UD0∂rDrD=εD=-q^D,UD0rD→∞=0,UD0tD=0=pmDtD=0=pvDtD=0=0.Equation (22) is the linearized line-sink model in triple-porosity reservoirs.


### 3.4.3. Laplace Transformation

Take the following Laplace transformation:(23)UD0¯=∫0+∞UD0e-utDdtD,where u is the Laplace transformation variable.Then (22) becomes(24)d2UD0¯drD2+1rDdUD0¯drD=ωfuUD0¯+ωmupmD¯+ωvupvD¯,(25)ωmupmD¯+λmfpmD¯-UD0¯=0,(26)ωvupvD¯+λvfpvD¯-UD0¯=0,(27)limεD→0⁡rDdUD0¯drDrD=εD=-q^D¯,(28)UD0¯rD→∞=0.


### 3.4.4. Line-Sink Solution

Equations (25) and (26) can be rewritten as(29)pmD¯=λmfωmu+λmfUD0¯,pvD¯=λvfωvu+λvfUD0¯.Substituting (29) into (24) yields(30)d2UD0¯drD2+1rDdUD0¯drD=ωf+λmfωmu+λmfωm+λvfωvu+λvfωvuUD0¯.If we define(31)gu=ωf+λmfωmu+λmfωm+λvfωvu+λvfωvu,then (30) can be expressed as follows:(32)d2UD0¯drD2+1rDdUD0¯drD=guUD0¯.The general solution of (32) can be easily obtained as follows:(33)UD0¯=c1K0gurD+c2I0gurD,where K0 and I0 are Bessel functions and c1 and c2 are unknown coefficients which can be determined by corresponding boundary conditions.If the line-sink is located at the origin of coordinates (center of the triple-porosity reservoir), then the dimensionless radial distancerD in (31) can be calculated by the following equation:(34a)rD=xD2+yD2.If the line-sink is located at(xw,yw) instead, then rD is calculated by the following equation:(34b)rD=xD-xwD2+yD-ywD2,where(35)xD=xL,(36)yD=yL,(37)xwD=xwL,(38)ywD=ywL.Substitution of (33) into the inner boundary condition (27) yields(39)limεD→0⁡gurD·-c1K1gurD+c2I1gurDrD=εD=-q^D¯.According to the properties of Bessel’s function, which arelimx→0⁡xK1(x)→1 and limx→0⁡xI1(x)→0, (39) becomes (40)c1=q^D¯.Similarly, with the outer boundary condition given in (28) and the properties of Bessel’s function, limx→∞⁡K0(x)→0 and limx→∞⁡I0(x)→∞, the coefficient in (33), c2, should satisfy the following equation:(41)c2=0.With (40) and (41), we can get the final form of (33) as follows:(42)UD0¯=q^D¯K0gurD.Equation (42) is the basic zero-order perturbation solution for a line-sink in a stress-sensitive triple-porosity reservoir. As mentioned above, zero-order perturbation solution is enough to approximate the exact solution of (15) to (20); that is, UD¯≈UD0¯. With the definition of q^D in Table 1, (42) can be expressed as(43)UD¯=q^¯qK0gurD.Equation (43) can be rewritten as(44)UD¯=q^¯qΩxD,yD;xwD,ywD,where(45)ΩxD,yD;xwD,ywD=K0gurD.Equation (44) is the basic line-sink solution for stress-sensitive triple-porosity reservoirs. With the basic line-sink solution and the superposition principle, we can obtain the pressure responses caused by multistage fractured horizontal wells in stress-sensitive triple-porosity reservoirs.


## 3.4.1. Pedrosa’s Linearization

The stress-sensitivity of natural fracture permeability makes the above seepage model strongly nonlinear which cannot be solved analytically. Here, following Pedrosa Jr. [25], we introduced the following equation to alleviate the nonlinearity of (8) and (11):(14)pfDrD,tD=-1γDln⁡1-γDUDrD,tD,where UD is an intermediate dimensionless parameter which is also a function of radial distance and time.With (14), (8)~(13) become(15)∂2UD∂rD2+1rD∂UD∂rD=11-γDUDωf∂UD∂tD+ωm∂pmD∂tD+ωv∂pvD∂tD,(16)ωm∂pmD∂tD+λmfpmD+1γDln⁡1-γDUD=0,(17)1-ωf-ωm∂pvD∂tD+λvfpvD+1γDln⁡1-γDUD=0,(18)limεD→0⁡rD∂UD∂rDrD=εD=-q^D,(19)UDrD→∞=0,(20)UDtD=0=pmDtD=0=pvDtD=0=0.


## 3.4.2. Perturbation Technique

According to the perturbation theory, the termsUD, -(1/γD)ln⁡[1-γDUD(rD,tD)], and 1/(1-γDUD(rD,tD)) in (15) through (20) can be expanded as power series in dimensionless permeability modulus, and we can get the following equations:(21)UD=UD0+γDUD1+γD2UD2+⋯,-1γDln⁡1-γDUDrD,tD=UDrD,tD+12γDUD2rD,tD+⋯,11-γDUDrD,tD=1+γDUDrD,tD+γD2UDrD,tD+⋯Given that the dimensionless permeability modulusγD is usually small (much smaller than 1), the zero-order approximate solution can satisfy the requirements of engineering precision, so (15)~(20) become(22)∂2UD0∂rD2+1rD∂UD0∂rD=ωf∂UD0∂tD+ωm∂pmD∂tD+ωv∂pvD∂tD,ωm∂pmD∂tD+λmfpmD-UD0=0,1-ωf-ωm∂pvD∂tD+λvfpvD-UD0=0,limεD→0⁡rD∂UD0∂rDrD=εD=-q^D,UD0rD→∞=0,UD0tD=0=pmDtD=0=pvDtD=0=0.Equation (22) is the linearized line-sink model in triple-porosity reservoirs.


## 3.4.3. Laplace Transformation

Take the following Laplace transformation:(23)UD0¯=∫0+∞UD0e-utDdtD,where u is the Laplace transformation variable.Then (22) becomes(24)d2UD0¯drD2+1rDdUD0¯drD=ωfuUD0¯+ωmupmD¯+ωvupvD¯,(25)ωmupmD¯+λmfpmD¯-UD0¯=0,(26)ωvupvD¯+λvfpvD¯-UD0¯=0,(27)limεD→0⁡rDdUD0¯drDrD=εD=-q^D¯,(28)UD0¯rD→∞=0.


## 3.4.4. Line-Sink Solution

Equations (25) and (26) can be rewritten as(29)pmD¯=λmfωmu+λmfUD0¯,pvD¯=λvfωvu+λvfUD0¯.Substituting (29) into (24) yields(30)d2UD0¯drD2+1rDdUD0¯drD=ωf+λmfωmu+λmfωm+λvfωvu+λvfωvuUD0¯.If we define(31)gu=ωf+λmfωmu+λmfωm+λvfωvu+λvfωvu,then (30) can be expressed as follows:(32)d2UD0¯drD2+1rDdUD0¯drD=guUD0¯.The general solution of (32) can be easily obtained as follows:(33)UD0¯=c1K0gurD+c2I0gurD,where K0 and I0 are Bessel functions and c1 and c2 are unknown coefficients which can be determined by corresponding boundary conditions.If the line-sink is located at the origin of coordinates (center of the triple-porosity reservoir), then the dimensionless radial distancerD in (31) can be calculated by the following equation:(34a)rD=xD2+yD2.If the line-sink is located at(xw,yw) instead, then rD is calculated by the following equation:(34b)rD=xD-xwD2+yD-ywD2,where(35)xD=xL,(36)yD=yL,(37)xwD=xwL,(38)ywD=ywL.Substitution of (33) into the inner boundary condition (27) yields(39)limεD→0⁡gurD·-c1K1gurD+c2I1gurDrD=εD=-q^D¯.According to the properties of Bessel’s function, which arelimx→0⁡xK1(x)→1 and limx→0⁡xI1(x)→0, (39) becomes (40)c1=q^D¯.Similarly, with the outer boundary condition given in (28) and the properties of Bessel’s function, limx→∞⁡K0(x)→0 and limx→∞⁡I0(x)→∞, the coefficient in (33), c2, should satisfy the following equation:(41)c2=0.With (40) and (41), we can get the final form of (33) as follows:(42)UD0¯=q^D¯K0gurD.Equation (42) is the basic zero-order perturbation solution for a line-sink in a stress-sensitive triple-porosity reservoir. As mentioned above, zero-order perturbation solution is enough to approximate the exact solution of (15) to (20); that is, UD¯≈UD0¯. With the definition of q^D in Table 1, (42) can be expressed as(43)UD¯=q^¯qK0gurD.Equation (43) can be rewritten as(44)UD¯=q^¯qΩxD,yD;xwD,ywD,where(45)ΩxD,yD;xwD,ywD=K0gurD.Equation (44) is the basic line-sink solution for stress-sensitive triple-porosity reservoirs. With the basic line-sink solution and the superposition principle, we can obtain the pressure responses caused by multistage fractured horizontal wells in stress-sensitive triple-porosity reservoirs.


## 4. Pressure Responses for Multistage Fractured Horizontal Wells in Stress-Sensitive Triple-Porosity Reservoirs

Figure4 shows the schematic of a horizontal well with multiple (a total of M) hydraulic fractures in a stress-sensitive triple-porosity reservoir. The y-axis is set to be aligned with the horizontal wellbore. The intersection of the ith (i=1,2,…,M) hydraulic fracture and the y-axis is represented by (0,yi), and the distance between every two adjacent fractures is Δyi.Figure 4
Schematic of discrete fracture segments.To successfully obtain the pressure responses caused by the multiple fractures with the line-sink solution derived in the above section, it is necessary to discretize the multiple hydraulic fractures. As shown in Figure4, each hydraulic fracture is divided into 2N segments along its length. The coordinate of the midpoint (i.e., discrete node) of the jth segment on the ith fracture is denoted by (Xi,j,Yi,j), and the coordinates of the two corresponding endpoints are denoted by (xi,j,yi,j) and (xi,j+1,yi,j+1).With the discretization scheme shown in Figure4, the coordinates of endpoints of discrete segments can be calculated by the following equations:(46)xi,j=-N-j+1NXfLi,yi,j=yi,1≤j≤N,xi,j=j-N-1NXfRi,yi,j=yi,N+1≤j≤2N+1.The coordinates of discrete nodes (midpoint of each segment) can be calculated by the following equations:(47)Xi,j=-2N-2j+12NXfLi,Yi,j=yi,1≤j≤N,Xi,j=2j-N-12NXfRi,Yi,j=yi,N+1≤j≤2N,where XfLi is the length of the left wing of ith fracture, m; XfRi is the length of the right wing of ith fracture, m. In the model presented in this paper, XfLi and XfRi can be different and can vary from fracture to fracture.The distance betweenyi and yi-1 is(48)Δyi=yi-yi-1.It is obvious that the flux strength is different along the fracture length. However, in the same discrete segment (assuming the discrete segment is small enough), the flux strength can be considered as constant. If we useqi,j to represent the flux density per unit length, then, with the superposition principle, the pressure response at (xD,yD) in the reservoir caused by the discrete segment (i,j) can be obtained by integrating the basic line-sink along the discrete segment, which is(49)UDi,j¯xD,yD=q-i,juq∫xi,jxi,j+1ΩxD,yD;xwD,ywDdxw.With the definitions ofxD and xwD, that is, (35) and (37), (49) can be changed into(50)UDi,j¯xD,yD=q-i,juLq∫xDi,jxDi,j+1ΩxD,yD;xwD,ywDdxwD.If we define dimensionless flux density per unit length as(51)qDi,jtD=qi,jtDLqand take Laplace transformation of (51),(52)q-Di,ju=q-i,juLq.Equation (50) can be rewritten as(53)UDi,j¯xD,yD=q-Di,j∫xDi,jxDi,j+1ΩxD,yD;xwD,ywDdxwD.The pressure response at(xD,yD) caused by 2N×M segments can be obtained by applying the principle of superposition over all discrete fracture segments:(54)UD¯xD,yD=∑i=1M∑j=12NUDi,j¯xD,yD.Thus, the pressure response at the discrete segment(χ,κ)(1≤χ≤M,1≤κ≤2N) can be obtained as follows:(55)UD¯xDχ,κ,yDχ,κ=∑i=1M∑j=12NUDi,j¯xDχ,κ,yDχ,κ.Since the permeability of hydraulic fractures is always much higher than the original reservoir permeability, the pressure drop in hydraulic fractures is much smaller than the pressure drop caused by oil flow in the reservoir. Thus, the hydraulic fractures can be considered as infinitely conductive, meaning the pressure in hydraulic fractures is equal to the bottom-hole pressure:(56)UD¯XDχ,κ,YDχ,κ=UwD¯.Combining (55) and (56), we can get(57)∑i=1M∑j=12NUDi,j¯XDχ,κ,YDχ,κ=UwD¯.The subscriptsχ and κ in (57) denote that this equation is written for discrete segment (χ,κ). By writing (57) for all discrete segments, that is, letting (χ=1,2,…,M,κ=1,2,…,2N), we can get 2N×M equations. It should be noted that there are 2N×M+1 unknowns in these 2N×M equations, which are UwD¯ and q-Di,j(i=1,2,…,M,j=1,2,…,2N). To compose a closed system of equations, the assumption of constant total flow rate can give the following equation: (58)∑i=1M∑j=12Nqi,jΔXi,j=q,where ΔXi,j is the length of the discrete segment (i,j) and can be calculated by(59)ΔXi,j=Xi,j+1-Xi,j.Taking Laplace transformation of (58), one can get(60)∑i=1M∑j=12Nq-i,jΔXi,j=qu.The dimensionless form of (60) is(61)∑i=1M∑j=12Nq-Di,jΔXDi,j=1u,where(62)ΔXDi,j=ΔXi,jL.Equations (57) and (61) represent a system of (2N×M+1) linear algebraic equations relating (2N×M+1) unknowns, which can be solved by direct methods, such as Gaussian elimination method, Gauss-Jordan reduction method.U
w
D
¯ calculated by (57) and (61) does not take into account the effects of wellbore storage and skin. According to van Everdingen [26] and Kucuk and Ayestaran [27], the effects of wellbore storage and skin can be incorporated in the calculated bottom-hole pressure response by(63)UwDS¯=uUwD¯+Su+CDu2uUwD¯+S,where CD is the dimensionless wellbore storage coefficient and S is the skin factor defined by the following, respectively: (64)CD=C2πϕfCft+ϕmCmt+ϕvCvthL2,(65)S=2πkfihqsfμΔps,where C is the wellbore storage coefficient, m3/Pa; qsf is the sandface flow rate of the multistage fractured horizontal well in the triple-porosity reservoir, m3/s; Δps is an extra pressure drop near the wellbore, Pa.U
w
D
S
¯ in (63) is the bottom-hole pressure response obtained in the Laplace domain, and by Stehfest [28] numerical inversion algorithm we can calculate the pressure responses UwDS(tD) in real time domain. Then, with (66), we can obtain the bottom-hole pressure response for multistage fractured horizontal wells in triple-porosity reservoirs incorporating the stress-sensitivity of natural fracture system:(66)pwDtD=-1γDln⁡1-γDUwDStD.


## 5. Results and Discussion

In this section, a horizontal well with three hydraulic fractures (i.e.,M=3) in a triple-porosity reservoir is considered, and the model introduced above is applied to investigate the dimensionless pressure responses of this well-reservoir configuration.


### 5.1. Flow Periods and Effect of Stress-Sensitivity

Figure5 presents the dimensionless pressure responses and corresponding pressure derivatives of a multistage fractured horizontal well in an infinitely large triple-porosity reservoir with consideration of stress-sensitivity of natural fractures. It can be observed from the figure that as many as ten flow periods can be identified at different time scales. Different flow periods and their corresponding characteristics are explained as follows.Figure 5
Log-log plots of dimensionless pressures and pressure derivatives of a multistage fractured horizontal well in an infinitely large triple-porosity reservoir.Period 1. It is the early-time wellbore storage period. In this period, both pressure curves and their derivative curves exhibit unit slope on the log-log plots.Period 2. It is the transition flow period after wellbore storage period. The pressure derivative curves exhibit a “hump,” representing the effect of skin.Period 3. It is the first linear flow period. This period corresponds to the linear flow in the direction perpendicular to hydraulic fracture length as illustrated in Figure 6(a). The pressure derivative curves exhibit a straight line with slope equal to “1/2.”Figure 6
(a) First linear flow period. (b) First pseudoradial flow period. (c) Second linear flow period. (d) Late-time compound pseudoradial flow period.
(a)
(b)
(c)
(d)Period 4. It is the first pseudoradial flow period. This period can occur when the distances between adjacent fractures are relatively large. Oil flows in the reservoir in a way as shown in Figure 6(b), and pseudoradial flow occurs around each fracture. This period corresponds to a nearly horizontal line with a value of “1/(2M)” in the derivative curve, where M is the number of hydraulic fractures.Period 5. It is the second linear flow period. This period corresponds to linear flow perpendicular to the horizontal wellbore as shown in Figure 6(c). Another straight line with a slope of about “1/2” can be observed in the derivative curve. When the permeability of natural fracture system remains constant during the whole production, corresponding to the no stress-sensitivity case, the slope of the derivative curves during this period is “1/2.” When taking into account the stress-sensitivity of natural fracture system, the corresponding pressure derivative curve slightly deviates from the one-half-slope line and exhibits a slight upward tendency.Period 6. It is the second pseudoradial flow period, corresponding to radial flow in natural fracture system. During this period, multiple hydraulic fractures and the horizontal wellbore behave like an enlarged wellbore. The derivative curves exhibit a horizontal line with an upward tendency depending on the value of dimensionless permeability modulus. In the case that γD=0, the derivative curves are horizontal with a value of 0.5 on the y-axis. With the increase of the value of γD, the derivative curves turn upward gradually.Period 7. It is the first interporosity flow period. This period reflects mass transfer between vugs and natural fracture system. Due to the oil flowing from vugs to natural fractures, the bottom-hole pressure drops slower and an obvious “dip” can be observed in the derivative curves.Period 8. It is the third pseudoradial flow period, corresponding to compound radial flow in natural fracture system and vugs. Due to the existence of stress-sensitivity of natural fracture system, the pressure derivative curves are no longer horizontal. The position of the derivative curves becomes higher with larger value of γD.Period 9. It is the second interporosity flow period. This period reflects oil flowing from matrix to natural fracture system. Another characteristic “dip” can be observed in the derivative curves.Period 10. It is the late-time compound pseudoradial flow period. This period reflects compound radial flow in natural fracture system, matrix, and vugs as shown in Figure 6(d). Similarly, the derivative curves exhibit upward tendencies due to the existence of stress-sensitivity of natural fracture system.


### 5.2. Effect of Skin Factor

Figure7 shows the effect of skin factor on dimensionless pressure curves and corresponding derivative curves. As stated above, the skin factor mainly affects the shape of type curves during the transition period (Period 2) after wellbore storage period. With the increase of the value of skin factor, the position of the “hump” in the derivative curves becomes higher, representing larger pressure drop in the formation.Figure 7
Type curves with different skin factors.


### 5.3. Effect of Storativity Ratio of Natural Fracture System

Figure8 shows the effect of storativity ratio of natural fracture system on dimensionless pressure and pressure derivative curves. It is obvious in Figure 8 that when the storativity ratio of matrix (i.e., ωm) and other parameters are kept constant, the storativity ratio of natural fracture system ωf mainly affects the transient pressure dynamics during the first linear flow period, the first pseudoradial flow period, the second linear flow period, second pseudoradial flow period, and the first interporosity flow period. On the one hand, the position of the dimensionless pressure curves becomes higher during the abovementioned flowing periods with the decrease of ωf. On the other hand, with the decrease of ωf, the position of the pressure derivative curves during the first linear flow period, the first pseudoradial flow, and the second linear flow period becomes higher, and the first “dip” in pressure derivative curves becomes wider and deeper.Figure 8
Type curves with different storativity ratios of natural fracture system.


### 5.4. Effect of Storativity Ratio of Matrix

Figure9 shows the effect of storativity ratio of matrix on dimensionless pressure and derivative curves. It can be observed that, with fixed ωf, the storativity ratio of matrix ωm has primary effect on pressure dynamics during two interporosity flowing periods. With the decrease of the value of ωm, the fist “dip” in the pressure derivative curves becomes wider and deeper, while the second “dip” in the pressure derivative curves exhibits the opposite trend.Figure 9
Type curves with different storativity ratios of matrix.


### 5.5. Effect of Interporosity Flow Coefficient between Vugs and Natural Fractures

Figure10 shows the effect of interporosity flow coefficient between vugs and natural fractures (λvf) on dimensionless pressure and pressure derivative curves. As expected, λvf mainly affects the occurrence time of the first interporosity flow period, that is, the appearance time of the first “dip” in pressure derivative curves. The smaller the value of λvf is, the later the first “dip” can be observed in the pressure derivative curves. When λvf decreases to a certain extent, the third pseudoradial flow period (Period 8 in Figure 5) may not be observed in the type curves. In addition, when λvf increases to a certain extent, the second pseudoradial flow period (Period 6 in Figure 5) and even the second linear flow period (Period 5 in Figure 5) may be masked by the following first interporosity flow period.Figure 10
Type curves with differentλvf.


### 5.6. Effect of Interporosity Flow Coefficient between Matrix and Natural Fractures

Figure11 shows the effect of interporosity flow coefficient between matrix and natural fractures (λmf) on dimensionless pressure and pressure derivative curves. It can be seen that λvf mainly affects the occurrence time of the second interporosity flow period, that is, the appearance time of the second “dip” in pressure derivative curves. The smaller the value of λmf is, the later the second “dip” can be observed in the pressure derivative curves. With the increase of λmf, the reflection of the third pseudoradial flow period on the derivative curves (Period 8 in Figure 5) may not be observed. When λmf increases to a certain degree, the second “dip” may merge with the first “dip,” meaning that simultaneous interporosity mass transfers happen between natural fractures and matrix as well as natural fractures and vugs.Figure 11
Type curves with differentλmf.


### 5.7. Effect of Length of Hydraulic Fractures

Figure12 presents the effect of half-length of hydraulic fractures on type curves. For the convenience of discussion, XfLi=XfRi=Xf is assumed here. It can be observed that the length of hydraulic fracture has primary effect on the first linear flow period and the subsequent first pseudoradial flow period. With the increase of fracture length, the first linear flow period lasts longer with lower position of pressure derivative curves, and the first pseudoradial flow period occurs later with shorter duration. When the half-length of hydraulic fractures continues to increase, the horizontal line in the pressure derivative curves reflecting the first pseudoradial flow period around each fracture will gradually disappear.Figure 12
Type curves with differentXf.


### 5.8. Effect of Hydraulic Fracture Spacing

Figure13 shows the effect of the distance between adjacent hydraulic fractures. The case discussed in Figure 13 assumes equal fracture spacing for the convenience of discussion. It can be observed that the distance between adjacent hydraulic fractures Δyi mainly affects the first pseudoradial flow period and the subsequent second linear flow period. The larger the distance between adjacent fractures, the longer the duration of the first pseudoradial flow period, and the later the occurrence of the second linear flow period. Another point should be addressed here is that when Δyi increases to a certain degree, the second pseudoradial flow period may not be observed in the type curves; that is, the first interporosity flow period appears after the second linear flow period.Figure 13
Type curves with differentΔyi.


### 5.9. Flux Distribution

Figure14 presents the dimensionless flux distribution along hydraulic fractures at a given time based on the calculation results obtained with the model proposed in this paper. In the case presented in Figure 14, three hydraulic fractures with equal half-length and fracture spacing are considered, and each hydraulic fracture is discretized into ten segments. It can be observed that, for the same hydraulic fracture, the flux distribution is symmetrical with respect to the horizontal wellbore, and the flux density at fracture tips is larger than that in the middle of the fracture. In addition, for different hydraulic fractures, it can be found that the flux contribution from hydraulic fracture in the middle (Fracture 3 in Figure 14) is always smaller than that from outer fractures (Fractures 1 and 2 in Figure 14).Figure 14
Dimensionless flux distribution along hydraulic fractures.


## 5.1. Flow Periods and Effect of Stress-Sensitivity

Figure5 presents the dimensionless pressure responses and corresponding pressure derivatives of a multistage fractured horizontal well in an infinitely large triple-porosity reservoir with consideration of stress-sensitivity of natural fractures. It can be observed from the figure that as many as ten flow periods can be identified at different time scales. Different flow periods and their corresponding characteristics are explained as follows.Figure 5
Log-log plots of dimensionless pressures and pressure derivatives of a multistage fractured horizontal well in an infinitely large triple-porosity reservoir.Period 1. It is the early-time wellbore storage period. In this period, both pressure curves and their derivative curves exhibit unit slope on the log-log plots.Period 2. It is the transition flow period after wellbore storage period. The pressure derivative curves exhibit a “hump,” representing the effect of skin.Period 3. It is the first linear flow period. This period corresponds to the linear flow in the direction perpendicular to hydraulic fracture length as illustrated in Figure 6(a). The pressure derivative curves exhibit a straight line with slope equal to “1/2.”Figure 6
(a) First linear flow period. (b) First pseudoradial flow period. (c) Second linear flow period. (d) Late-time compound pseudoradial flow period.
(a)
(b)
(c)
(d)Period 4. It is the first pseudoradial flow period. This period can occur when the distances between adjacent fractures are relatively large. Oil flows in the reservoir in a way as shown in Figure 6(b), and pseudoradial flow occurs around each fracture. This period corresponds to a nearly horizontal line with a value of “1/(2M)” in the derivative curve, where M is the number of hydraulic fractures.Period 5. It is the second linear flow period. This period corresponds to linear flow perpendicular to the horizontal wellbore as shown in Figure 6(c). Another straight line with a slope of about “1/2” can be observed in the derivative curve. When the permeability of natural fracture system remains constant during the whole production, corresponding to the no stress-sensitivity case, the slope of the derivative curves during this period is “1/2.” When taking into account the stress-sensitivity of natural fracture system, the corresponding pressure derivative curve slightly deviates from the one-half-slope line and exhibits a slight upward tendency.Period 6. It is the second pseudoradial flow period, corresponding to radial flow in natural fracture system. During this period, multiple hydraulic fractures and the horizontal wellbore behave like an enlarged wellbore. The derivative curves exhibit a horizontal line with an upward tendency depending on the value of dimensionless permeability modulus. In the case that γD=0, the derivative curves are horizontal with a value of 0.5 on the y-axis. With the increase of the value of γD, the derivative curves turn upward gradually.Period 7. It is the first interporosity flow period. This period reflects mass transfer between vugs and natural fracture system. Due to the oil flowing from vugs to natural fractures, the bottom-hole pressure drops slower and an obvious “dip” can be observed in the derivative curves.Period 8. It is the third pseudoradial flow period, corresponding to compound radial flow in natural fracture system and vugs. Due to the existence of stress-sensitivity of natural fracture system, the pressure derivative curves are no longer horizontal. The position of the derivative curves becomes higher with larger value of γD.Period 9. It is the second interporosity flow period. This period reflects oil flowing from matrix to natural fracture system. Another characteristic “dip” can be observed in the derivative curves.Period 10. It is the late-time compound pseudoradial flow period. This period reflects compound radial flow in natural fracture system, matrix, and vugs as shown in Figure 6(d). Similarly, the derivative curves exhibit upward tendencies due to the existence of stress-sensitivity of natural fracture system.


## 5.2. Effect of Skin Factor

Figure7 shows the effect of skin factor on dimensionless pressure curves and corresponding derivative curves. As stated above, the skin factor mainly affects the shape of type curves during the transition period (Period 2) after wellbore storage period. With the increase of the value of skin factor, the position of the “hump” in the derivative curves becomes higher, representing larger pressure drop in the formation.Figure 7
Type curves with different skin factors.


## 5.3. Effect of Storativity Ratio of Natural Fracture System

Figure8 shows the effect of storativity ratio of natural fracture system on dimensionless pressure and pressure derivative curves. It is obvious in Figure 8 that when the storativity ratio of matrix (i.e., ωm) and other parameters are kept constant, the storativity ratio of natural fracture system ωf mainly affects the transient pressure dynamics during the first linear flow period, the first pseudoradial flow period, the second linear flow period, second pseudoradial flow period, and the first interporosity flow period. On the one hand, the position of the dimensionless pressure curves becomes higher during the abovementioned flowing periods with the decrease of ωf. On the other hand, with the decrease of ωf, the position of the pressure derivative curves during the first linear flow period, the first pseudoradial flow, and the second linear flow period becomes higher, and the first “dip” in pressure derivative curves becomes wider and deeper.Figure 8
Type curves with different storativity ratios of natural fracture system.


## 5.4. Effect of Storativity Ratio of Matrix

Figure9 shows the effect of storativity ratio of matrix on dimensionless pressure and derivative curves. It can be observed that, with fixed ωf, the storativity ratio of matrix ωm has primary effect on pressure dynamics during two interporosity flowing periods. With the decrease of the value of ωm, the fist “dip” in the pressure derivative curves becomes wider and deeper, while the second “dip” in the pressure derivative curves exhibits the opposite trend.Figure 9
Type curves with different storativity ratios of matrix.


## 5.5. Effect of Interporosity Flow Coefficient between Vugs and Natural Fractures

Figure10 shows the effect of interporosity flow coefficient between vugs and natural fractures (λvf) on dimensionless pressure and pressure derivative curves. As expected, λvf mainly affects the occurrence time of the first interporosity flow period, that is, the appearance time of the first “dip” in pressure derivative curves. The smaller the value of λvf is, the later the first “dip” can be observed in the pressure derivative curves. When λvf decreases to a certain extent, the third pseudoradial flow period (Period 8 in Figure 5) may not be observed in the type curves. In addition, when λvf increases to a certain extent, the second pseudoradial flow period (Period 6 in Figure 5) and even the second linear flow period (Period 5 in Figure 5) may be masked by the following first interporosity flow period.Figure 10
Type curves with differentλvf.


## 5.6. Effect of Interporosity Flow Coefficient between Matrix and Natural Fractures

Figure11 shows the effect of interporosity flow coefficient between matrix and natural fractures (λmf) on dimensionless pressure and pressure derivative curves. It can be seen that λvf mainly affects the occurrence time of the second interporosity flow period, that is, the appearance time of the second “dip” in pressure derivative curves. The smaller the value of λmf is, the later the second “dip” can be observed in the pressure derivative curves. With the increase of λmf, the reflection of the third pseudoradial flow period on the derivative curves (Period 8 in Figure 5) may not be observed. When λmf increases to a certain degree, the second “dip” may merge with the first “dip,” meaning that simultaneous interporosity mass transfers happen between natural fractures and matrix as well as natural fractures and vugs.Figure 11
Type curves with differentλmf.


## 5.7. Effect of Length of Hydraulic Fractures

Figure12 presents the effect of half-length of hydraulic fractures on type curves. For the convenience of discussion, XfLi=XfRi=Xf is assumed here. It can be observed that the length of hydraulic fracture has primary effect on the first linear flow period and the subsequent first pseudoradial flow period. With the increase of fracture length, the first linear flow period lasts longer with lower position of pressure derivative curves, and the first pseudoradial flow period occurs later with shorter duration. When the half-length of hydraulic fractures continues to increase, the horizontal line in the pressure derivative curves reflecting the first pseudoradial flow period around each fracture will gradually disappear.Figure 12
Type curves with differentXf.


## 5.8. Effect of Hydraulic Fracture Spacing

Figure13 shows the effect of the distance between adjacent hydraulic fractures. The case discussed in Figure 13 assumes equal fracture spacing for the convenience of discussion. It can be observed that the distance between adjacent hydraulic fractures Δyi mainly affects the first pseudoradial flow period and the subsequent second linear flow period. The larger the distance between adjacent fractures, the longer the duration of the first pseudoradial flow period, and the later the occurrence of the second linear flow period. Another point should be addressed here is that when Δyi increases to a certain degree, the second pseudoradial flow period may not be observed in the type curves; that is, the first interporosity flow period appears after the second linear flow period.Figure 13
Type curves with differentΔyi.


## 5.9. Flux Distribution

Figure14 presents the dimensionless flux distribution along hydraulic fractures at a given time based on the calculation results obtained with the model proposed in this paper. In the case presented in Figure 14, three hydraulic fractures with equal half-length and fracture spacing are considered, and each hydraulic fracture is discretized into ten segments. It can be observed that, for the same hydraulic fracture, the flux distribution is symmetrical with respect to the horizontal wellbore, and the flux density at fracture tips is larger than that in the middle of the fracture. In addition, for different hydraulic fractures, it can be found that the flux contribution from hydraulic fracture in the middle (Fracture 3 in Figure 14) is always smaller than that from outer fractures (Fractures 1 and 2 in Figure 14).Figure 14
Dimensionless flux distribution along hydraulic fractures.


## 6. Conclusion

Based on the mathematical model and discussion in this paper, the following conclusions can be warranted:(1)
A mathematical model is presented for investigating transient pressure dynamics as well as flux distribution of multistage fractured horizontal wells in stress-sensitive triple-porosity reservoirs. The model can better describe fluid flow in naturally fractured reservoirs with multiple types of storage spaces and complex well type.(2)
Semianalytical solution is developed in the Laplace domain by the integral transformation method and discretization of hydraulic fractures. With the semianalytical solution and the numerical inversion algorithm proposed by Stehfest, dimensionless pressure responses and corresponding pressure derivatives as well as flux distribution among all the fractures can be obtained in real time domain. Computation results prove that the model presented in this paper is stable and accurate.(3)
Analysis of transient pressure and derivative responses indicates that as many as ten flow periods can be identified for a multistage fractured horizontal well in triple-porosity reservoirs. Major factors affecting the transient pressure responses include skin factor, storativity ratios of different porous media type, interporosity flow coefficients, and parameters related to the geometry and placement of hydraulic fractures. With different combinations of these parameters, the reflection of some flow periods may be masked in the pressure derivative curves.(4)
Two characteristic “dips” can be observed in the pressure derivative curves for triple-porosity reservoirs, reflecting the mass transfer between natural fracture-vugs and natural fracture-matrix, separately.(5)
The stress-sensitivity of natural fracture system results in larger pressure drop during intermediate and late flowing periods which is reflected by upward tendencies in both dimensionless pressure curves and corresponding derivative curves.(6)
Analysis of flux distribution among multiple hydraulic fractures indicates that the flux contribution from outer fractures is higher, and for a given fracture the flux contribution from fracture tips is higher.(7)
The presented model can be used to interpret pressure signal for fractured horizontal wells in triple-porosity reservoirs and provide some important dynamic parameters for reservoir development.



---

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