# Pressure Transient Analysis for a Horizontal Well in Heterogeneous Carbonate Reservoirs Using a Linear Composite Model

Authors: Yong-Gang Duan; Ke-Yi Ren; Quan-Tang Fang; Ming-Qiang Wei; Morteza Dejam; Wei-Hua Chen
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Carbonate reservoirs usually have strong anisotropy. Oil and gas recovery from fractured reservoirs is highly challenging due to complicated mechanisms involved in production from these reservoirs. A horizontal well completed in these reservoirs may extend through multiple zones, including homogeneous, dual-porosity, and triple-porosity formations. Traditional well test models assume that the entire length of a horizontal or multilateral well remains in the same formation with uniform properties. A well test model for pressure transient analysis of horizontal wells extending through a carbonate reservoir consisting of natural fractures, rock matrix, and vugs with different properties is presented in this study. The focus of this study is on dual-porosity (fracture-matrix) and triple-porosity (fracture-matrix-vug) reservoirs, considering the pseudosteady interporosity flows from rock matrix and vugs into fractures. A multizone triple-porosity model was established and solved by using the point source function, Green’s function, and coupling of multiple reservoir sections. The corresponding type curves were developed, and sensitivity analysis was carried out. The type curves of flow stage division reveal that a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix)/triple-porosity (fracture-vug-matrix), and homogeneous sections identifies the stages of pseudosteady interporosity flow from matrix and vug into fracture, fracture pseudoradial flow, system linear flow, system pseudoradial flow, and pseudosteady flow occur in sequence. The greater the difference of permeability between the dual-porosity/triple-porosity section and the two homogeneous sections, the more obvious the interporosity flow on the pressure derivative curve. This approach satisfies the need for pressure transient analysis for a horizontal well that traverses two or more regions with distinct properties in heterogeneous carbonate reservoirs.


---

## Body

## 1. Introduction

Carbonate reservoirs play a very important role in the world’s gas and oil production. Most of the world’s giant fields produce from naturally fractured and vuggy carbonate reservoirs that have complex pore structure. Drilling horizontal wells would be challenging because of unfavorable geomechanical properties [1, 2]. The majority of carbonate reservoirs are very heterogeneous consisting of different pore classes as well as being fractured. One of the most prominent problems of carbonate reservoirs is the strong heterogeneity and multiporosity network characteristic, which makes the fluid flow system complicated and causes relatively low oil recovery. Many mature theories and practices of clastic rocks in reservoir development cannot be applied to carbonate reservoirs. It is generally considered that carbonate reservoirs have their own characteristics, that is, well-developed connected fracture systems, most of which also contain tiny pore systems connected with fracture systems. Unlike the conventional single-porosity sandstone reservoirs, these kinds of reservoirs need to consider the permeability and porosity of the rock mass as well as the exchange of quality between the fractures and matrix in addition to considering the permeability and porosity of the fracture network as the main seepage channel. Different from the common conceptual dual-porosity coal reservoirs where gas is mostly stored in the coal matrix and Darcy fluid flow occurs in the natural fracture system [3, 4] and shale gas reservoirs where gas is trapped in poorly connected micropores by absorption on and in particles of organic materials and clay minerals in the matrix of the host rock [5], fractures have high permeability and low porosity, while matrix has low permeability and high porosity. The low flow in the matrix mainly comes from the fluid exchange between matrix and fracture network. The presence of fractures and vugs significantly influences the dynamic of fluid flow in carbonate reservoirs and the transient pressure analysis of wells drilled in these reservoirs.Barenblatt et al. [6] pointed out that it is not promising to study the nonsteady-state seepage in fractured rocks by assuming a system of fractures, and therefore, they proposed the concept of dual porosity for fluid flow in naturally fractured reservoirs for the first time. They mentioned that the interconnected fractures and matrix rock in fractured reservoirs as two continuous media are spatially superimposed on each other. The permeability of matrix rock is far less than that of the fracture system, but the porosity of matrix rock is far higher than that of the fracture system. Based on this hypothesis, Warren and Root [7] developed an idealized model for studying the characteristic behavior of a porous medium, which contains regions that contribute significantly to the pore volume of the system but contribute negligibly to the flow capacity. Motivated by the early model by Warren and Root, the unsteady flow in dual-porosity media was mathematically modeled and solved analytically [8, 9]. Then, with the development of modern well testing methods, the effect of matrix-block-size distribution on pressure transient response for naturally fractured reservoirs was discussed [10–13]. Considering both pseudosteady state and unsteady-state fluid transfer from matrix blocks into fracture network, Raghavan and Ohaeri [14] delineated the characteristics of a well flowing at a constant pressure in a naturally fractured reservoir.However, most of the fractured carbonate reservoirs in the world have more complex pore structure. Vugs contribute to both effective porosity and permeability. Matrix and vugs do not have the same interaction with the fracture network [15]. Abdassah and Ershaghi [16] first proposed a triple-porosity/single-permeability model, which considers an unsteady-state interporosity flow between the fractures and two types of matrix blocks and a primary flow only through the fracture network. Jalali and Ershaghi [17] analyzed the unsteady pressure behavior of such reservoirs. New dual-fracture models [18] considering two types of flow units were proposed to differentiate between the microfractures and the macrofractures. Then, a triple porosity model [15] was built to address high secondary porosity and mainly vuggy porosity in naturally fractured reservoirs.Horizontal and multifractured horizontal wells can improve the flow conductivity and enhance well production and recovery factor by increasing the exposure of wells to formation [19, 20]. Traditional pressure transient models assume that the entire length of a horizontal or multilateral well remains in the same formation with uniform properties. In reality, most horizontal and multilateral wells extend through the reservoirs consisting of layers, sections, or pockets with different properties [21]. Bowman and Crawford [22] presented a sufficient number of tables to calculate the pressure transient in a linear semi-infinite water-driven reservoir, where the properties in the oil and water zones were different. Bixel et al. [23] proposed a detailed treatment of the pressure transient behavior for a well located near a linear discontinuity. Strltsova and Mckinley [24] discussed the buildup pattern for linear discontinuities and the effect of flow time on quantitative assessment of reservoir parameters. Ambastha et al. [25] employed analytical solutions to obtain pressure-transient behavior for a line-source well in a composite reservoir. Based on Ambastha’s model, an analytical solution for a three-zone laterally composite reservoir model [26] was found in the Laplace–Fourier space, and Kuchuk and Tarek [27] extended this model to multizone case, explaining the changes in porosity, permeability, and compressibility. Based on the material balance equation and considering the nonlinearity caused by gas desorption and diffusion, and the complex flow of MFHW, Wei et al. [28, 29] developed a numerical method for multifractured horizontal wells in shale gas reservoirs. Medeiros et al. [21] presented a semianalytical model for the pressure-transient analysis of horizontal wells in linear composite reservoirs and provided the detailed solutions for a two-zone, laterally composite reservoir. If the reservoir is divided into more than two sections, their solution methodology will not work. The reason is that the coupling of middle sections of multizone reservoir is much more complicated than two-zone reservoir. Dejam et al. [30–32] used semianalytical solutions for pressure transient analysis of a hydraulically fractured vertical well and multifractured horizontal well in a bounded dual-porosity reservoir. As carbonate reservoirs usually have strong anisotropy, a horizontal well completed in these reservoirs may extend through multiple zones, including homogeneous, dual-porosity, and triple-porosity formations. Traditional well test models are inapplicable here. Therefore, it is necessary to build new models to satisfy the need for well testing in a horizontal well in heterogeneous carbonate reservoirs.This study is developed based on Medeiros et al.’s model and works for a horizontal well extending through a multizone reservoir. Considering different combinations of natural fractures, rock matrix, and vugs in carbonate reservoirs, various physical and mathematical models are established and solved by using the point source function, Green’s function, and coupling of multiple reservoir sections. The corresponding type curves are obtained by computer programming, and sensitivity analysis are carried out.


## 2. Solution Methodology

First, in order to solve the problem under study, where a horizontal well traverses two or more regions with distinct properties, Green’s function formulation of the pressure transient solution for a locally homogeneous reservoir section with arbitrary flux conditions at the boundaries [21] is applied. Then, the coupling of different sections of the reservoir is explained.


### 2.1. Source Functions for the Inner and Outer Boundaries

Green’s function for a rectangular parallelepiped satisfies the adjoint diffusion equation with its flux vanishing at the boundaries of the domain. This Green’s function corresponds to the point source solution in the Laplace transform domain for a rectangular parallelepiped one in Figure1 [33].Figure 1
Rectangular parallelepiped representing the reservoir block for the point source solution [34].(1)
G
¯
M
D
,
M
D
′
,
s
=
μ
2
K
l
x
e
D
h
D
c
h
u
y
˜
D
1
+
c
h
u
y
˜
D
2
u
s
h
u
y
˜
e
D
+
2
∑
k
=
1
∞
cos
k
π
x
D
x
e
D
cos
k
π
x
D
′
x
e
D
c
h
ε
k
y
˜
D
1
+
c
h
ε
k
y
˜
D
2
ε
k
s
h
u
y
˜
e
D
+
2
∑
n
=
1
∞
cos
n
π
z
D
z
e
D
cos
n
π
z
D
′
z
e
D
c
h
ε
n
y
˜
D
1
+
c
h
ε
n
y
˜
D
2
ε
n
s
h
u
y
˜
e
D
+
2
∑
k
=
1
∞
cos
k
π
x
D
x
e
D
cos
k
π
x
D
′
x
e
D
c
h
ε
k
,
n
y
˜
D
1
+
c
h
ε
k
,
n
y
˜
D
2
ε
k
,
n
s
h
u
y
˜
e
D
,Where(2)hD=zeD=hlKKz,(3)y˜D1=yeD−yD−yD′,(4)y˜D2=yeD−yD+yD′,(5)εn=u+nπhD2,(6)εk=u+kπxeD2,(7)εk,n=u+kπxeD2+nπhD2,where D is the domain consisting of porous medium, eD is the outer boundary of the domain, l is the characteristic reference length in the system, K is the uniform permeability, Kz is the principal permeability in the z direction, and μ is the viscosity of a slightly compressible fluid of constant compressibility.In equation (1) and equations (5)–(7), the new parameter u was introduced to extend the solution into the dual-porosity or triple-porosity idealization of naturally fractured reservoirs [33, 34] and is given by(8)u=s,for homogeneous reservoirs,sfs,for naturally fractured reservoirs.Appropriate expressions forfs are related to the type of formation [34, 35]. This paper mainly considers the conventional dual-porosity and the triple-porosity models, where the seepage type is triporate-uniphase parallel flow. The corresponding expressions for fs are given, respectively, as follows:(9)fs=sw1−w+λs1−w+λ,(10)fs=wNf+wvλvswv+λv+wmλmswm+λm,where w is the dimensionless storativity and λ is the dimensionless transfer coefficient. The subscript Nf denotes the natural fracture system, v refers to the vug, m represents the matrix, and s stands for the Laplace variable.Based on Medeiros and Ozkan’s model, the inner boundary,Bw, is a line equal to the length of the horizontal well, Lh. Then, the source function in the Laplace domain, S¯α, given by Medeiros and Ozkan can be expressed as(11)S¯wi=∫xwi−Lhi/2xwi+Lhi/2G¯xD,yD,zD,xD′,ywiD,zwiD,sdx′,where xwi, ywi, and zwi are the coordinates of the midpoint and Lhi is the length of the ith segment of the horizontal well. Substituting G¯ given by equation (1) into equation (11) and evaluating the integral yield,(12)S¯wiMD,s=CLhichuy˜D1i+chuy˜D2iushuy˜eD+2C∑k=1∞chεky˜D1i+chεky˜D2iεkshεkuy˜eD∫xwi−Lhi/2xwi+Lhi/2coskπxxecoskπx′xedx′+2CLhi∑n=1∞cosnπzDzeDcosnπzwiDzeDchεny˜D1i+chεny˜D2iεnshεnuy˜eD+4CLhi∑n=1∞cosnπzDzeDcosnπzwiDzeD×∑j=1∞chεj,ny˜D1i+chεj,ny˜D2iεj,nshεj,nuy˜eD∫xwi−Lhi/2xwi+Lhi/2cosjπxxecosjπx′xedx′,where(13)C=μ2KlxeDhD.The outer boundary consists of the planar surfaces of the rectangular parallelepiped representing the domain. In this study, the source function for a segment on the boundary perpendicular to thex-axis (a planar source in the y-z plane) is given by(14)S¯ej=∫yej−ypj/2yej+ypj/2∫zej−zpj/2zej+zpj/2G¯xD,yD,zD,xwD,y′,z′,sdz′dy′,where yej and zej are the coordinates of the midpoint and ypj and zpj are the lengths of the sides of the planar segment.Substituting equation (1) into equation (14) and evaluating the internal result in(15)S¯ejMD,s=Czpj∫yej−ypj/2yej+ypj/2chuy˜D1j+chuy˜D2jushuy˜eDdy′+2Czpj∑k=1∞coskπxDxeDcoskπxwDxeD∫yej−ypj/2yej+ypj/2chεky˜D1j+chεky˜D2jεkshεky˜eDdy′+2C∑n=1∞∫zej−zpj/2zej+zpj/2cosnπzzecosnπz′zedz′×∫yej−ypj/2yej+ypj/2chεny˜D1j+chεny˜D2jεnshεny˜eDdy′+4C∫zej−zpj/2zej+zpj/2∑n=1∞cosnπzzecosnπz′zedz′×∑k=1∞coskπxDxeDcoskπxwDxeD∫yej−ypj/2yej+ypj/2chεk,ny˜D1j+chεk,ny˜D2jεk,nshεk,ny˜eDdy′.


### 2.2. Coupling of Multiple Reservoir Sections

Medeiros et al. [21] developed a semianalytical model for the pressure-transient analysis of horizontal wells in linear composite reservoirs and provided the detailed solution for a two-zone reservoir. In their study, all section boundaries except the interface between sections are assumed impermeable. Physically, this case corresponds to a horizontal well that penetrates a closed and rectangular reservoir with two homogeneous sections of different properties. However, their model cannot provide solution if the horizontal well traverses more than two sections. The reason is that the coupling of middle sections of a multizone reservoir is much more complicated than a two-zone reservoir.To demonstrate the coupling of multiple reservoir sections by using the solution given in equation (1), two rectangular reservoir sections in series in the x-direction are considered as the basis for the derivation. Then, the models for three reservoir sections and even multiple reservoir sections models are developed. For the sake of simplicity, the discretization scheme shown in Figures 2 and 3, which includes three horizontal well segments and two interface segments between the sections, is used.Figure 2
Three-section reservoir with a horizontal well.Figure 3
Discretization scheme to demonstrate the coupling of three reservoir sections penetrated by a horizontal well.From equation (1), the pressure drop at the center of a well or interface segment i in section k is given by(16)Δp¯MDi,s=∑j=12q˜¯wjksS¯iwjkMDi,s+∑l=12q˜¯elksS¯ielkM,Dis,where MDi indicates the midpoint of segment i, q˜¯wj and q˜¯el are the fluxes on all segments of the inner and outer boundaries, respectively. Using equation (16) at the center of all well and interface segments shown in Figure 3 yields the following set of 14 linear equations with 28 unknowns:(17)q˜¯w1kS¯wiw1k+q˜¯w2kS¯wiw2k+q˜¯e1kS¯wie1k+q˜¯e2kS¯wie2k−Δp¯wik=0,for i = 1, 2, k = 1, 3,(18)q˜¯w1kS¯wiw1k+q˜¯w2kS¯wiw2k+q˜¯e1k1S¯wie1k1+q˜¯e2k1S¯wie2k1+q˜¯e1k2S¯wie1k2+q˜¯e2k2S¯wie2k2−Δp¯wik=0,for i = 1, 2, k = 2,(19)q˜¯w1kS¯eiw1k+q˜¯w2kS¯eiw2k+q˜¯e1kS¯eie1k+q˜¯e2kS¯eie2k−Δp¯eik=0,for i = 1, 2, k = 1, 3,(20)q˜¯w12S¯eiw12k+q˜¯w22S¯eiw22k+q˜¯e121S¯eie12k1+q˜¯e221S¯eie22k1−q˜¯e122S¯eie12k2+q˜¯e222S¯eie22k2−Δp¯ei2k=0,for i, k = 1, 2.To equate the number of equations with the number of unknowns, the condition of continuity of pressure and flux at the interface between sections 1, 2, and 3 is implemented. The continuity of pressure and flux requires(21)Δp¯ei1=Δp¯ei21,Δp¯ei22=Δp¯ei3,fori=1,2,(22)q˜¯ei1=−q˜¯ei21,q˜¯ei22=−q˜¯ei3,fori=1,2.This reduces the number of unknowns by twelve. If it is also assumed that there is infinite conductivity in the wellbore, then five additional equations are obtained:(23)Δp¯w1k−Δp¯w2k=0,fork=1,2,3,(24)Δp¯w21−Δp¯w12=0,Δp¯w22−Δp¯w13=0.Finally, the mass balance requires that the sum of fluxes entering the wellbore is equal to the production rate at the sandfaceq:(25)∑k=13∑j=12q˜¯wjkLhj=qs.The linear system defined by equation (17)–(25) now has 16 equations with 16 unknowns.For models addressing multiple reservoir sections, the discretization scheme shown in Figures4 and 5, which includes n horizontal well segments and n − 1 interface segments between the sections, is used.Figure 4
n-section reservoir with a horizontal well.Figure 5
Discretization scheme to demonstrate the coupling ofn reservoir sections penetrated by a horizontal well.Using equation (16) at the center of all well and interface segments shown in Figure 5 yields the following set of 2(3n − 2) linear equations with 4(3n − 2) unknowns:(26)q˜¯w1kS¯wiw1k+q˜¯w2kS¯wiw2k+q˜¯e1kS¯wie1k+q˜¯e2kS¯wie2k−Δp¯wik=0,for i = 1, 2, k = 1, n,(27)q˜¯w1kS¯wiw1k+q˜¯w2kS¯wiw2k+q˜¯e1k1S¯wie1k1+q˜¯e2k1S¯wie2k1+q˜¯e1k2S¯wie1k2+q˜¯e2k2S¯wie2k2−Δp¯wik=0,for i = 1, 2, k = 2, 3, …, n − 1,(28)q˜¯w1kS¯eiw1k+q˜¯w2kS¯eiw2k+q˜¯e1kS¯eie1k+q˜¯e2kS¯eie2k−Δp¯eik=0,for i = 1, 2, k = 1, n,(29)q˜¯w1kS¯eiw1km+q˜¯w2kS¯eiw2km+q˜¯e1k1S¯eie1km1+q˜¯e2k1S¯eie2km1−q˜¯e1k2S¯eie1km2+q˜¯e2k2S¯eie2km2−Δp¯eikm=0,for i, m = 1, 2, k = 2, 3, …, n − 1.To equate the number of equations with the number of unknowns, the condition of continuity of pressure and flux at the interface between sections 1, 2, …,n is used. The continuity of pressure and flux requires(30)Δp¯ei1=Δp¯ei21,Δp¯eik2=Δp¯eik+11,Δp¯ein−12=Δp¯ein,fori=1,2,k=2,3,…,n−1,(31)q˜¯ei1=−q˜¯ei21,q˜¯eik2=−q˜¯eik+11,q˜¯ein−12=−q˜¯ein,fori=1,2,k=2,3,…,n−1.This reduces the number of unknowns by 4n. If it is also assumed that there is infinite conductivity in the wellbore, then 2n − 1 additional equations are obtained:(32)Δp¯w1k−Δp¯w2k=0,fork=1,2,3,…,n,(33)Δp¯w2k−Δp¯w1k+1=0,fork=1,2,3,…,n−1.Finally, mass balance requires that the sum of fluxes entering the wellbore is equal to the production rate at the sandfaceq:(34)∑k=1n∑j=12q˜¯wjkLhj=qs.The linear system defined by equations (26)–(34) now has 2(3n−1) equations with 2(3n−1) unknowns.


### 2.3. Solution of the System of Linear Equations

The system of linear equations, given by equations (16)–(25), is in the matrix-vector form, Ax = b. The components of the coefficient matrix, A = {A1, A2}, are given by(35)A2=−1000000−1000000000000000000−1000000−1000000000000000000−1000000−11−1000001−1000001−1000001−1000001−1000000,A1=S¯w1w11S¯w1w21S¯w1e11S¯w1e21000000S¯w2w11S¯w2w21S¯w2e11S¯w2e21000000S¯e1w11S¯e1w21−S¯e1e11+S¯e1e1211−S¯e1e21+S¯e1e2211−S¯e1w121−S¯e1w221S¯e1e1212S¯e1e221200S¯e2w11S¯e2w21−S¯e2e11+S¯e2e1211−S¯e2e21+S¯e2e2211−S¯e2w121−S¯e2w221S¯e2e1212S¯e2e22120000S¯w1e121S¯w1e221S¯w1w12S¯w1w22−S¯w1e122−S¯w1e2220000S¯w2e121S¯w2e221S¯w2w12S¯w2w22−S¯w2e122−S¯w2e2220000S¯e1e1221S¯e1e2221S¯e1w122S¯e1w222−S¯e1e1222+S¯e1e13−S¯e1e2222+S¯e1e23S¯e1w13S¯e1w2300S¯e2e1221S¯e2e2221S¯e2w122S¯e1w222−S¯e2e1222+S¯e2e13−S¯e2e2222+S¯e2e23S¯e2w13S¯e2w23000000S¯w1e13S¯w1e23S¯w1w13S¯w1w23000000S¯w2e13S¯w2e23S¯w2w13S¯w2w2300000000000000000000000000000000000000000000000000Lh11Lh2100Lh12Lh2200Lh13Lh23.The solution vector,x = [x1, x2], has the following components:(36)x1=q˜¯w11q˜¯w21q˜¯e121q˜¯e221q˜¯w12q˜¯w12q˜¯e13q˜¯e23q˜¯w13q˜¯w13,x2=Δp¯w11Δp¯w21Δp¯w12Δp¯w22Δp¯w13Δp¯w23.The components of the right-hand-side vector,b = [b1, b2], are(37)b1=0000000000,b2=00000q/s.The system of linear equations, given by equations (26)–(34), is in the matrix-vector form, Ax = b. The components of the coefficient matrix, A = {A1, A2}, are given by(38)A2=−1−100−1−100⋱4n−3rows−1−11−11−11−1⋱⋱2n−3rows1−11−11−10⋯2n−2columns⋯0,A1=S¯w1w1k−1S¯w1w2k−1S¯w1e1k−12S¯w1e2k−12S¯w2w1k−1S¯w2w2k−1S¯w2e1k−12S¯w2e2k−12S¯e1w1k−1S¯e1w2k−12−S¯e1e1k−122+S¯e1e1k11−S¯e1e2k−122+S¯e1e2k11−S¯e1w1k1−S¯e1w2k1S¯e1e1k12S¯e1e2k12S¯e2w1k−1S¯e2w2k−12−S¯e2e1k−122+S¯e2e1k11−S¯e2e2k−122+S¯e2e2k11−S¯e2w1k1−S¯e2w2k1S¯e2e1k12S¯e2e2k12S¯w1e1k1S¯w1e2k1S¯w1w1kS¯w1w2k−S¯w1e1k2−S¯w1e2k2S¯w2e1k1S¯w2e2k1S¯w2w1kS¯w2w2k−S¯w2e1k2−S¯w2e2k2⋯S¯e1e1k21S¯e1e2k21S¯e1w1k2S¯e1w2k2−S¯e1e1k22+S¯e1e1k+111−S¯e1e2k22+S¯e1e2k+111S¯e1w1k+11S¯e1w2k+11⋯S¯e2e1k21S¯e2e2k21S¯e2w1k2S¯e2w2k2−S¯e2e1k22+S¯e2e1k+111−S¯e2e2k22+S¯e2e2k+111S¯e2w1k+11S¯e2w2k+11S¯w1e1k+11S¯w1e2k+11S¯w1w1k+1S¯w1w2k+1S¯w2e1k+11S¯w2e2k+11S¯w2w1k+1S¯w2w2k+10⋱0⋯Lh1k−1Lh2k−100Lh1kLh2k00Lh1k+1Lh2k+1⋯.The solution vector,x = [x1, x2], has the following components:(39)x1=q˜¯w11q˜¯w21⋯q˜¯e1k1q˜¯e2k1q˜¯w1kq˜¯w1k⋯q˜¯e1nq˜¯e2nq˜¯w1nq˜¯w1n,fork=2,3,…,n−1,x2=Δp¯w11Δp¯w21⋯Δp¯w1kΔp¯w2k⋯Δp¯w1nΔp¯w2n,fork=2,3,…,n−1.The components of the right-hand-side vector,b = [b1, b2], are(40)b1=00⋯00⋯001×4n−2,b2=00⋯00q/s1×2n.The solution of this system of linear equations can be accomplished by using any matrix inverter. It is, however, in the Laplace transform domain and needs to be numerically inverted into the real-time domain. As it is common in the petroleum engineering literature, the Stehfest algorithm [36] is used for the numerical inversion of Laplace transforms in this study.


## 2.1. Source Functions for the Inner and Outer Boundaries

Green’s function for a rectangular parallelepiped satisfies the adjoint diffusion equation with its flux vanishing at the boundaries of the domain. This Green’s function corresponds to the point source solution in the Laplace transform domain for a rectangular parallelepiped one in Figure1 [33].Figure 1
Rectangular parallelepiped representing the reservoir block for the point source solution [34].(1)
G
¯
M
D
,
M
D
′
,
s
=
μ
2
K
l
x
e
D
h
D
c
h
u
y
˜
D
1
+
c
h
u
y
˜
D
2
u
s
h
u
y
˜
e
D
+
2
∑
k
=
1
∞
cos
k
π
x
D
x
e
D
cos
k
π
x
D
′
x
e
D
c
h
ε
k
y
˜
D
1
+
c
h
ε
k
y
˜
D
2
ε
k
s
h
u
y
˜
e
D
+
2
∑
n
=
1
∞
cos
n
π
z
D
z
e
D
cos
n
π
z
D
′
z
e
D
c
h
ε
n
y
˜
D
1
+
c
h
ε
n
y
˜
D
2
ε
n
s
h
u
y
˜
e
D
+
2
∑
k
=
1
∞
cos
k
π
x
D
x
e
D
cos
k
π
x
D
′
x
e
D
c
h
ε
k
,
n
y
˜
D
1
+
c
h
ε
k
,
n
y
˜
D
2
ε
k
,
n
s
h
u
y
˜
e
D
,Where(2)hD=zeD=hlKKz,(3)y˜D1=yeD−yD−yD′,(4)y˜D2=yeD−yD+yD′,(5)εn=u+nπhD2,(6)εk=u+kπxeD2,(7)εk,n=u+kπxeD2+nπhD2,where D is the domain consisting of porous medium, eD is the outer boundary of the domain, l is the characteristic reference length in the system, K is the uniform permeability, Kz is the principal permeability in the z direction, and μ is the viscosity of a slightly compressible fluid of constant compressibility.In equation (1) and equations (5)–(7), the new parameter u was introduced to extend the solution into the dual-porosity or triple-porosity idealization of naturally fractured reservoirs [33, 34] and is given by(8)u=s,for homogeneous reservoirs,sfs,for naturally fractured reservoirs.Appropriate expressions forfs are related to the type of formation [34, 35]. This paper mainly considers the conventional dual-porosity and the triple-porosity models, where the seepage type is triporate-uniphase parallel flow. The corresponding expressions for fs are given, respectively, as follows:(9)fs=sw1−w+λs1−w+λ,(10)fs=wNf+wvλvswv+λv+wmλmswm+λm,where w is the dimensionless storativity and λ is the dimensionless transfer coefficient. The subscript Nf denotes the natural fracture system, v refers to the vug, m represents the matrix, and s stands for the Laplace variable.Based on Medeiros and Ozkan’s model, the inner boundary,Bw, is a line equal to the length of the horizontal well, Lh. Then, the source function in the Laplace domain, S¯α, given by Medeiros and Ozkan can be expressed as(11)S¯wi=∫xwi−Lhi/2xwi+Lhi/2G¯xD,yD,zD,xD′,ywiD,zwiD,sdx′,where xwi, ywi, and zwi are the coordinates of the midpoint and Lhi is the length of the ith segment of the horizontal well. Substituting G¯ given by equation (1) into equation (11) and evaluating the integral yield,(12)S¯wiMD,s=CLhichuy˜D1i+chuy˜D2iushuy˜eD+2C∑k=1∞chεky˜D1i+chεky˜D2iεkshεkuy˜eD∫xwi−Lhi/2xwi+Lhi/2coskπxxecoskπx′xedx′+2CLhi∑n=1∞cosnπzDzeDcosnπzwiDzeDchεny˜D1i+chεny˜D2iεnshεnuy˜eD+4CLhi∑n=1∞cosnπzDzeDcosnπzwiDzeD×∑j=1∞chεj,ny˜D1i+chεj,ny˜D2iεj,nshεj,nuy˜eD∫xwi−Lhi/2xwi+Lhi/2cosjπxxecosjπx′xedx′,where(13)C=μ2KlxeDhD.The outer boundary consists of the planar surfaces of the rectangular parallelepiped representing the domain. In this study, the source function for a segment on the boundary perpendicular to thex-axis (a planar source in the y-z plane) is given by(14)S¯ej=∫yej−ypj/2yej+ypj/2∫zej−zpj/2zej+zpj/2G¯xD,yD,zD,xwD,y′,z′,sdz′dy′,where yej and zej are the coordinates of the midpoint and ypj and zpj are the lengths of the sides of the planar segment.Substituting equation (1) into equation (14) and evaluating the internal result in(15)S¯ejMD,s=Czpj∫yej−ypj/2yej+ypj/2chuy˜D1j+chuy˜D2jushuy˜eDdy′+2Czpj∑k=1∞coskπxDxeDcoskπxwDxeD∫yej−ypj/2yej+ypj/2chεky˜D1j+chεky˜D2jεkshεky˜eDdy′+2C∑n=1∞∫zej−zpj/2zej+zpj/2cosnπzzecosnπz′zedz′×∫yej−ypj/2yej+ypj/2chεny˜D1j+chεny˜D2jεnshεny˜eDdy′+4C∫zej−zpj/2zej+zpj/2∑n=1∞cosnπzzecosnπz′zedz′×∑k=1∞coskπxDxeDcoskπxwDxeD∫yej−ypj/2yej+ypj/2chεk,ny˜D1j+chεk,ny˜D2jεk,nshεk,ny˜eDdy′.


## 2.2. Coupling of Multiple Reservoir Sections

Medeiros et al. [21] developed a semianalytical model for the pressure-transient analysis of horizontal wells in linear composite reservoirs and provided the detailed solution for a two-zone reservoir. In their study, all section boundaries except the interface between sections are assumed impermeable. Physically, this case corresponds to a horizontal well that penetrates a closed and rectangular reservoir with two homogeneous sections of different properties. However, their model cannot provide solution if the horizontal well traverses more than two sections. The reason is that the coupling of middle sections of a multizone reservoir is much more complicated than a two-zone reservoir.To demonstrate the coupling of multiple reservoir sections by using the solution given in equation (1), two rectangular reservoir sections in series in the x-direction are considered as the basis for the derivation. Then, the models for three reservoir sections and even multiple reservoir sections models are developed. For the sake of simplicity, the discretization scheme shown in Figures 2 and 3, which includes three horizontal well segments and two interface segments between the sections, is used.Figure 2
Three-section reservoir with a horizontal well.Figure 3
Discretization scheme to demonstrate the coupling of three reservoir sections penetrated by a horizontal well.From equation (1), the pressure drop at the center of a well or interface segment i in section k is given by(16)Δp¯MDi,s=∑j=12q˜¯wjksS¯iwjkMDi,s+∑l=12q˜¯elksS¯ielkM,Dis,where MDi indicates the midpoint of segment i, q˜¯wj and q˜¯el are the fluxes on all segments of the inner and outer boundaries, respectively. Using equation (16) at the center of all well and interface segments shown in Figure 3 yields the following set of 14 linear equations with 28 unknowns:(17)q˜¯w1kS¯wiw1k+q˜¯w2kS¯wiw2k+q˜¯e1kS¯wie1k+q˜¯e2kS¯wie2k−Δp¯wik=0,for i = 1, 2, k = 1, 3,(18)q˜¯w1kS¯wiw1k+q˜¯w2kS¯wiw2k+q˜¯e1k1S¯wie1k1+q˜¯e2k1S¯wie2k1+q˜¯e1k2S¯wie1k2+q˜¯e2k2S¯wie2k2−Δp¯wik=0,for i = 1, 2, k = 2,(19)q˜¯w1kS¯eiw1k+q˜¯w2kS¯eiw2k+q˜¯e1kS¯eie1k+q˜¯e2kS¯eie2k−Δp¯eik=0,for i = 1, 2, k = 1, 3,(20)q˜¯w12S¯eiw12k+q˜¯w22S¯eiw22k+q˜¯e121S¯eie12k1+q˜¯e221S¯eie22k1−q˜¯e122S¯eie12k2+q˜¯e222S¯eie22k2−Δp¯ei2k=0,for i, k = 1, 2.To equate the number of equations with the number of unknowns, the condition of continuity of pressure and flux at the interface between sections 1, 2, and 3 is implemented. The continuity of pressure and flux requires(21)Δp¯ei1=Δp¯ei21,Δp¯ei22=Δp¯ei3,fori=1,2,(22)q˜¯ei1=−q˜¯ei21,q˜¯ei22=−q˜¯ei3,fori=1,2.This reduces the number of unknowns by twelve. If it is also assumed that there is infinite conductivity in the wellbore, then five additional equations are obtained:(23)Δp¯w1k−Δp¯w2k=0,fork=1,2,3,(24)Δp¯w21−Δp¯w12=0,Δp¯w22−Δp¯w13=0.Finally, the mass balance requires that the sum of fluxes entering the wellbore is equal to the production rate at the sandfaceq:(25)∑k=13∑j=12q˜¯wjkLhj=qs.The linear system defined by equation (17)–(25) now has 16 equations with 16 unknowns.For models addressing multiple reservoir sections, the discretization scheme shown in Figures4 and 5, which includes n horizontal well segments and n − 1 interface segments between the sections, is used.Figure 4
n-section reservoir with a horizontal well.Figure 5
Discretization scheme to demonstrate the coupling ofn reservoir sections penetrated by a horizontal well.Using equation (16) at the center of all well and interface segments shown in Figure 5 yields the following set of 2(3n − 2) linear equations with 4(3n − 2) unknowns:(26)q˜¯w1kS¯wiw1k+q˜¯w2kS¯wiw2k+q˜¯e1kS¯wie1k+q˜¯e2kS¯wie2k−Δp¯wik=0,for i = 1, 2, k = 1, n,(27)q˜¯w1kS¯wiw1k+q˜¯w2kS¯wiw2k+q˜¯e1k1S¯wie1k1+q˜¯e2k1S¯wie2k1+q˜¯e1k2S¯wie1k2+q˜¯e2k2S¯wie2k2−Δp¯wik=0,for i = 1, 2, k = 2, 3, …, n − 1,(28)q˜¯w1kS¯eiw1k+q˜¯w2kS¯eiw2k+q˜¯e1kS¯eie1k+q˜¯e2kS¯eie2k−Δp¯eik=0,for i = 1, 2, k = 1, n,(29)q˜¯w1kS¯eiw1km+q˜¯w2kS¯eiw2km+q˜¯e1k1S¯eie1km1+q˜¯e2k1S¯eie2km1−q˜¯e1k2S¯eie1km2+q˜¯e2k2S¯eie2km2−Δp¯eikm=0,for i, m = 1, 2, k = 2, 3, …, n − 1.To equate the number of equations with the number of unknowns, the condition of continuity of pressure and flux at the interface between sections 1, 2, …,n is used. The continuity of pressure and flux requires(30)Δp¯ei1=Δp¯ei21,Δp¯eik2=Δp¯eik+11,Δp¯ein−12=Δp¯ein,fori=1,2,k=2,3,…,n−1,(31)q˜¯ei1=−q˜¯ei21,q˜¯eik2=−q˜¯eik+11,q˜¯ein−12=−q˜¯ein,fori=1,2,k=2,3,…,n−1.This reduces the number of unknowns by 4n. If it is also assumed that there is infinite conductivity in the wellbore, then 2n − 1 additional equations are obtained:(32)Δp¯w1k−Δp¯w2k=0,fork=1,2,3,…,n,(33)Δp¯w2k−Δp¯w1k+1=0,fork=1,2,3,…,n−1.Finally, mass balance requires that the sum of fluxes entering the wellbore is equal to the production rate at the sandfaceq:(34)∑k=1n∑j=12q˜¯wjkLhj=qs.The linear system defined by equations (26)–(34) now has 2(3n−1) equations with 2(3n−1) unknowns.


## 2.3. Solution of the System of Linear Equations

The system of linear equations, given by equations (16)–(25), is in the matrix-vector form, Ax = b. The components of the coefficient matrix, A = {A1, A2}, are given by(35)A2=−1000000−1000000000000000000−1000000−1000000000000000000−1000000−11−1000001−1000001−1000001−1000001−1000000,A1=S¯w1w11S¯w1w21S¯w1e11S¯w1e21000000S¯w2w11S¯w2w21S¯w2e11S¯w2e21000000S¯e1w11S¯e1w21−S¯e1e11+S¯e1e1211−S¯e1e21+S¯e1e2211−S¯e1w121−S¯e1w221S¯e1e1212S¯e1e221200S¯e2w11S¯e2w21−S¯e2e11+S¯e2e1211−S¯e2e21+S¯e2e2211−S¯e2w121−S¯e2w221S¯e2e1212S¯e2e22120000S¯w1e121S¯w1e221S¯w1w12S¯w1w22−S¯w1e122−S¯w1e2220000S¯w2e121S¯w2e221S¯w2w12S¯w2w22−S¯w2e122−S¯w2e2220000S¯e1e1221S¯e1e2221S¯e1w122S¯e1w222−S¯e1e1222+S¯e1e13−S¯e1e2222+S¯e1e23S¯e1w13S¯e1w2300S¯e2e1221S¯e2e2221S¯e2w122S¯e1w222−S¯e2e1222+S¯e2e13−S¯e2e2222+S¯e2e23S¯e2w13S¯e2w23000000S¯w1e13S¯w1e23S¯w1w13S¯w1w23000000S¯w2e13S¯w2e23S¯w2w13S¯w2w2300000000000000000000000000000000000000000000000000Lh11Lh2100Lh12Lh2200Lh13Lh23.The solution vector,x = [x1, x2], has the following components:(36)x1=q˜¯w11q˜¯w21q˜¯e121q˜¯e221q˜¯w12q˜¯w12q˜¯e13q˜¯e23q˜¯w13q˜¯w13,x2=Δp¯w11Δp¯w21Δp¯w12Δp¯w22Δp¯w13Δp¯w23.The components of the right-hand-side vector,b = [b1, b2], are(37)b1=0000000000,b2=00000q/s.The system of linear equations, given by equations (26)–(34), is in the matrix-vector form, Ax = b. The components of the coefficient matrix, A = {A1, A2}, are given by(38)A2=−1−100−1−100⋱4n−3rows−1−11−11−11−1⋱⋱2n−3rows1−11−11−10⋯2n−2columns⋯0,A1=S¯w1w1k−1S¯w1w2k−1S¯w1e1k−12S¯w1e2k−12S¯w2w1k−1S¯w2w2k−1S¯w2e1k−12S¯w2e2k−12S¯e1w1k−1S¯e1w2k−12−S¯e1e1k−122+S¯e1e1k11−S¯e1e2k−122+S¯e1e2k11−S¯e1w1k1−S¯e1w2k1S¯e1e1k12S¯e1e2k12S¯e2w1k−1S¯e2w2k−12−S¯e2e1k−122+S¯e2e1k11−S¯e2e2k−122+S¯e2e2k11−S¯e2w1k1−S¯e2w2k1S¯e2e1k12S¯e2e2k12S¯w1e1k1S¯w1e2k1S¯w1w1kS¯w1w2k−S¯w1e1k2−S¯w1e2k2S¯w2e1k1S¯w2e2k1S¯w2w1kS¯w2w2k−S¯w2e1k2−S¯w2e2k2⋯S¯e1e1k21S¯e1e2k21S¯e1w1k2S¯e1w2k2−S¯e1e1k22+S¯e1e1k+111−S¯e1e2k22+S¯e1e2k+111S¯e1w1k+11S¯e1w2k+11⋯S¯e2e1k21S¯e2e2k21S¯e2w1k2S¯e2w2k2−S¯e2e1k22+S¯e2e1k+111−S¯e2e2k22+S¯e2e2k+111S¯e2w1k+11S¯e2w2k+11S¯w1e1k+11S¯w1e2k+11S¯w1w1k+1S¯w1w2k+1S¯w2e1k+11S¯w2e2k+11S¯w2w1k+1S¯w2w2k+10⋱0⋯Lh1k−1Lh2k−100Lh1kLh2k00Lh1k+1Lh2k+1⋯.The solution vector,x = [x1, x2], has the following components:(39)x1=q˜¯w11q˜¯w21⋯q˜¯e1k1q˜¯e2k1q˜¯w1kq˜¯w1k⋯q˜¯e1nq˜¯e2nq˜¯w1nq˜¯w1n,fork=2,3,…,n−1,x2=Δp¯w11Δp¯w21⋯Δp¯w1kΔp¯w2k⋯Δp¯w1nΔp¯w2n,fork=2,3,…,n−1.The components of the right-hand-side vector,b = [b1, b2], are(40)b1=00⋯00⋯001×4n−2,b2=00⋯00q/s1×2n.The solution of this system of linear equations can be accomplished by using any matrix inverter. It is, however, in the Laplace transform domain and needs to be numerically inverted into the real-time domain. As it is common in the petroleum engineering literature, the Stehfest algorithm [36] is used for the numerical inversion of Laplace transforms in this study.


## 3. Results and Discussion

### 3.1. Homogeneous-Dual-Porosity (Fracture-Matrix)-Homogeneous Model

It is assumed that a horizontal well traverses a reservoir consisting of homogeneous, dual-porosity (fracture-matrix), and homogeneous sections with the same permeability. The physical model is illustrated in Figure6.Figure 6
Homogeneous-dual-porosity (fracture-matrix)-homogeneous model.In Figure7, the solid curves represent the type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections. The triangles represent the type curves for a dual-porosity infinite-conductivity model. The stages of pseudosteady interporosity flow from matrix into fracture, fracture pseudoradial flow, system linear flow, system pseudoradial flow, and pseudosteady flow occur in sequence. In the early stage of pseudosteady interporosity flow from matrix into fracture, a valley appears on the pressure derivative curve. Due to the effect of the homogeneous sections, the valley of the three-section model is shallower. Then, the fracture pseudoradial flow happens with a constant pressure derivative value. In the stage of system linear flow, the derivative curve exhibits a straight line with a slope of 0.5, and thereafter, the system pseudoradial flow occurs. Finally, there is a pseudosteady flow stage with a slope of 1 on the pressure derivative curve.Figure 7
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections (solid curves). The triangles represent the type curves for a dual-porosity infinite conductivity model.It can be seen from Figure8 that, in the early stage, due to the interporosity flow from matrix into fracture, the flow rate of the two segments of wellbore in the dual-porosity section is larger than the flow rate of the two segments of wellbore in the homogeneous sections. With the advent of pseudoradial flow regime, the rates for all six segments of the wellbore become equal. At late times, the flux distribution tends to be stable and the segment closer to the heel of the horizontal well gets higher flux.Figure 8
Flow rate per segment for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections.Based on this model, keeping the other parameters unchanged, when the permeability of two homogeneous sections decreases to 10 mD, the pressure drop and pressure derivative curves become different, as shown in Figure9.Figure 9
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections (solid curves). The triangles represent the type curves for the model with the decreased permeability for the two homogeneous sections.In Figure9, the solid curves represent the original type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections. The triangles represent the type curves for the model with the decreased permeability for the two homogeneous sections. When the permeability of two homogeneous sections decreases from 100 mD to 10 mD, the values of pressure drop and pressure derivative increase. In addition, the stage of the interporosity flow from matrix into fracture becomes easy to recognize in the early period. The valley becomes deeper. Other flow regimes remain unchanged.In Figure10, the flow rate of the dual-porosity section is larger than that of the two homogeneous sections because the permeability of section 2 is higher than that of sections 1 and 3. In the early stage, due to the influence of interporosity flow, the wellbore flow in the second section gradually decreases until the stage of the fracture pseudoradial flow appears. With the increase of time, the rates for all six segments of the wellbore reach a steady value.Figure 10
Flux distribution along the horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections (k1 = k3 = 10 mD).Based on the model, when the permeability of sections 1 and 2 keeps unchanged and the permeability of the homogeneous reservoir in section 3 increases from 10 mD to 40 mD and then to 100 mD, the wellbore pressure drops and pressure derivative values decreases as shown in Figure11. When the permeability of homogeneous section 3 becomes closer to that of section 2 (dual-porosity), the dip of the interporosity flow stage becomes shallower.Figure 11
Type curves of flow stage division when the permeability of section 3 increases from 10 mD to 40 mD and then to 100 mD.Figure12 shows the flux distribution along the horizontal well when the permeability of section 3 is 40 mD. The flow rate of the two segments of wellbore in section 2 with the highest permeability is greater than that of the two segments of wellbore in other two sections.Figure 12
Flow rate per segment when the permeability of section 3 is 40 mD.


### 3.2. Homogeneous-Triple-Porosity (Fracture-Vug-Matrix)-Homogeneous Model

It is assumed that a horizontal well traverses a reservoir consisting of homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections with the same permeability. The physical model is illustrated in Figure13.Figure 13
Homogeneous-triple-porosity (fracture-vug-matrix)-homogeneous model.In Figure14, the solid curves represent the type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections. The triangles represent the type curves for a triple-porosity infinite conductivity model. The stages of pseudosteady interporosity flows from vug and matrix into fracture, fracture pseudoradial flow, system linear flow, system pseudoradial flow, and pseudosteady flow occur sequentially. In the early stage of pseudosteady interporosity flow from vug and matrix into fracture, two valleys appear on the pressure derivative curve. Due to the effect of the homogeneous sections, the valleys of the three-section model are shallower. Then, the fracture pseudoradial flow happens with a constant pressure derivative value. In the stage of system linear flow, the derivative curve exhibits a straight line with a slope of 0.5, and thereafter, the system pseudoradial flow occurs. Finally, there is a pseudosteady flow stage with a slope of 1 on the pressure derivative curve.Figure 14
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections (solid curves). The triangles represent the type curves for a triple-porosity infinite conductivity model.It can be observed from the flux distribution in Figure15 that due to the two interporosity flows of the triple-porosity section in the early stage, the flow rate of the two segments of wellbore in the triple-porosity section is larger than the flow rate of the two segments of wellbore in the homogeneous sections. With the appearance of pseudoradial flow regime, the rates for all six segments of the wellbore become equal. At late times, the flux distribution tends to be stable, and the segment closer to the heel of the horizontal well gets higher flux.Figure 15
Flow rate per segment for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections.Based on this model, keeping the other parameters unchanged, when the permeability of two homogeneous sections decrease to 10 mD, the pressure drops and pressure derivative curves become different, as shown in Figure16.Figure 16
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections (solid curves). The triangles represent the type curves for the model with the decreased permeability for the two homogeneous sections.In Figure16, the solid curves represent the original type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections. The triangles represent the type curves for the model with the decreased permeability for the two homogeneous sections. When the permeability of homogeneous sections decreases from 100 mD to 10 mD, the values of pressure drop and pressure derivative increase and the stage of the interporosity flow from both vug and matrix into fracture becomes easier to recognize in the early period. The dips become deeper. Other flow regimes remain unchanged.It can be seen from Figure17 that the flow rate of the two segments of wellbore in section 2 is higher than that of the two segments of wellbore in other two sections. That is because the permeability of the triple-porosity section is higher than that of the first and third sections. In the early stage, due to the influence of interporosity flows, the wellbore flow in the second section gradually decreases until the stage of the fracture pseudoradial flow appears. With the increase of time, the rates for all six segments of the wellbore reach a steady value.Figure 17
Flux distribution along the horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections (k1 = k3 = 10 mD).


### 3.3. Homogeneous-Dual-Porosity (Fracture-Matrix)-Triple-Porosity (Fracture-Vug-Matrix) Model

It is assumed that a horizontal well traverses a reservoir consisting of homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug- matrix) sections with the same permeability. The physical model is illustrated in Figure18.Figure 18
Homogeneous-dual-porosity (fracture-matrix)-triple-porosity (fracture-vug-matrix) model.In Figure19, the solid curves represent a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug-matrix) sections. The triangles represent the type curves for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections. When a horizontal well traverses a reservoir consisting of homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug- matrix) sections with equal permeability, due to the influence of the first homogeneous section and the middle dual-porosity section, the interporosity flow from vug into fracture becomes weaker than the interporosity flow from matrix into fracture. The dip of the interporosity flow in dual-porosity section becomes easier to recognize; however, the general flow regimes remain unchanged.Figure 19
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug-matrix) sections (solid curves). The triangles represent the type curves for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections.In the flux distribution in Figure20, due to the two interporosity flows of the triple-porosity section and the interporosity flow of the dual-porosity section in the early stage, the flow rate of the two segments of wellbore in the triple-porosity section and the flow rate of the two segments of wellbore in the dual-porosity section are larger than the flow rate of the two segments of wellbore in the homogeneous section. With the appearance of pseudoradial flow regime, the rates for all six segments of the wellbore become equal. At late time, the flux distribution tends to be stable and the segment closer to the heel of the horizontal well gets higher flux.Figure 20
Flow rate per segment for a horizontal well traverses homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug-matrix) sections.Based on this model, when the permeability of homogeneous section decreases from 100 mD to 40 mD and then to 10 mD, the values of pressure drop and pressure derivative increase (Figure21). In addition, the stage of the interporosity flow from both vug and matrix into fracture system becomes harder to recognize at early times. The dip of the interporosity flow in dual-porosity section becomes shallower. Other flow regimes remain unchanged.Figure 21
Type curves of flow stage division when the permeability of section 1 increases from 10 mD to 40 mD and then to 100 mD.When the permeability of dual-porosity section decreases from 200 mD to 100 mD and then to 50 mD, the values of pressure drop and pressure derivative increase as it can be observed in Figure22. In addition, the stage of the interporosity flow from both vug and matrix into fracture becomes harder to recognize at early times. However, the dip of the interporosity flow in dual-porosity section becomes easier to identify. Other flow regimes remain unchanged.Figure 22
Type curves of flow stage division when the permeability of section 2 increases from 50 mD to 100 mD and then to 200 mD.Figure23 shows that when the permeability of triple-porosity section decreases from 200 mD to 100 mD and then to 50 mD, the values of pressure drop and pressure derivative increase. In addition, the stage of the interporosity flow from both vug and matrix into fracture becomes easier to recognize ate early times. Moreover, the dips of the interporosity flows in triple-porosity section become deeper. Other flow regimes remain unchanged.Figure 23
Type curves of flow stage division when the permeability of section 3 increases from 50 mD to 100 mD and then to 200 mD.


## 3.1. Homogeneous-Dual-Porosity (Fracture-Matrix)-Homogeneous Model

It is assumed that a horizontal well traverses a reservoir consisting of homogeneous, dual-porosity (fracture-matrix), and homogeneous sections with the same permeability. The physical model is illustrated in Figure6.Figure 6
Homogeneous-dual-porosity (fracture-matrix)-homogeneous model.In Figure7, the solid curves represent the type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections. The triangles represent the type curves for a dual-porosity infinite-conductivity model. The stages of pseudosteady interporosity flow from matrix into fracture, fracture pseudoradial flow, system linear flow, system pseudoradial flow, and pseudosteady flow occur in sequence. In the early stage of pseudosteady interporosity flow from matrix into fracture, a valley appears on the pressure derivative curve. Due to the effect of the homogeneous sections, the valley of the three-section model is shallower. Then, the fracture pseudoradial flow happens with a constant pressure derivative value. In the stage of system linear flow, the derivative curve exhibits a straight line with a slope of 0.5, and thereafter, the system pseudoradial flow occurs. Finally, there is a pseudosteady flow stage with a slope of 1 on the pressure derivative curve.Figure 7
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections (solid curves). The triangles represent the type curves for a dual-porosity infinite conductivity model.It can be seen from Figure8 that, in the early stage, due to the interporosity flow from matrix into fracture, the flow rate of the two segments of wellbore in the dual-porosity section is larger than the flow rate of the two segments of wellbore in the homogeneous sections. With the advent of pseudoradial flow regime, the rates for all six segments of the wellbore become equal. At late times, the flux distribution tends to be stable and the segment closer to the heel of the horizontal well gets higher flux.Figure 8
Flow rate per segment for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections.Based on this model, keeping the other parameters unchanged, when the permeability of two homogeneous sections decreases to 10 mD, the pressure drop and pressure derivative curves become different, as shown in Figure9.Figure 9
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections (solid curves). The triangles represent the type curves for the model with the decreased permeability for the two homogeneous sections.In Figure9, the solid curves represent the original type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections. The triangles represent the type curves for the model with the decreased permeability for the two homogeneous sections. When the permeability of two homogeneous sections decreases from 100 mD to 10 mD, the values of pressure drop and pressure derivative increase. In addition, the stage of the interporosity flow from matrix into fracture becomes easy to recognize in the early period. The valley becomes deeper. Other flow regimes remain unchanged.In Figure10, the flow rate of the dual-porosity section is larger than that of the two homogeneous sections because the permeability of section 2 is higher than that of sections 1 and 3. In the early stage, due to the influence of interporosity flow, the wellbore flow in the second section gradually decreases until the stage of the fracture pseudoradial flow appears. With the increase of time, the rates for all six segments of the wellbore reach a steady value.Figure 10
Flux distribution along the horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and homogeneous sections (k1 = k3 = 10 mD).Based on the model, when the permeability of sections 1 and 2 keeps unchanged and the permeability of the homogeneous reservoir in section 3 increases from 10 mD to 40 mD and then to 100 mD, the wellbore pressure drops and pressure derivative values decreases as shown in Figure11. When the permeability of homogeneous section 3 becomes closer to that of section 2 (dual-porosity), the dip of the interporosity flow stage becomes shallower.Figure 11
Type curves of flow stage division when the permeability of section 3 increases from 10 mD to 40 mD and then to 100 mD.Figure12 shows the flux distribution along the horizontal well when the permeability of section 3 is 40 mD. The flow rate of the two segments of wellbore in section 2 with the highest permeability is greater than that of the two segments of wellbore in other two sections.Figure 12
Flow rate per segment when the permeability of section 3 is 40 mD.


## 3.2. Homogeneous-Triple-Porosity (Fracture-Vug-Matrix)-Homogeneous Model

It is assumed that a horizontal well traverses a reservoir consisting of homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections with the same permeability. The physical model is illustrated in Figure13.Figure 13
Homogeneous-triple-porosity (fracture-vug-matrix)-homogeneous model.In Figure14, the solid curves represent the type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections. The triangles represent the type curves for a triple-porosity infinite conductivity model. The stages of pseudosteady interporosity flows from vug and matrix into fracture, fracture pseudoradial flow, system linear flow, system pseudoradial flow, and pseudosteady flow occur sequentially. In the early stage of pseudosteady interporosity flow from vug and matrix into fracture, two valleys appear on the pressure derivative curve. Due to the effect of the homogeneous sections, the valleys of the three-section model are shallower. Then, the fracture pseudoradial flow happens with a constant pressure derivative value. In the stage of system linear flow, the derivative curve exhibits a straight line with a slope of 0.5, and thereafter, the system pseudoradial flow occurs. Finally, there is a pseudosteady flow stage with a slope of 1 on the pressure derivative curve.Figure 14
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections (solid curves). The triangles represent the type curves for a triple-porosity infinite conductivity model.It can be observed from the flux distribution in Figure15 that due to the two interporosity flows of the triple-porosity section in the early stage, the flow rate of the two segments of wellbore in the triple-porosity section is larger than the flow rate of the two segments of wellbore in the homogeneous sections. With the appearance of pseudoradial flow regime, the rates for all six segments of the wellbore become equal. At late times, the flux distribution tends to be stable, and the segment closer to the heel of the horizontal well gets higher flux.Figure 15
Flow rate per segment for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections.Based on this model, keeping the other parameters unchanged, when the permeability of two homogeneous sections decrease to 10 mD, the pressure drops and pressure derivative curves become different, as shown in Figure16.Figure 16
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections (solid curves). The triangles represent the type curves for the model with the decreased permeability for the two homogeneous sections.In Figure16, the solid curves represent the original type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections. The triangles represent the type curves for the model with the decreased permeability for the two homogeneous sections. When the permeability of homogeneous sections decreases from 100 mD to 10 mD, the values of pressure drop and pressure derivative increase and the stage of the interporosity flow from both vug and matrix into fracture becomes easier to recognize in the early period. The dips become deeper. Other flow regimes remain unchanged.It can be seen from Figure17 that the flow rate of the two segments of wellbore in section 2 is higher than that of the two segments of wellbore in other two sections. That is because the permeability of the triple-porosity section is higher than that of the first and third sections. In the early stage, due to the influence of interporosity flows, the wellbore flow in the second section gradually decreases until the stage of the fracture pseudoradial flow appears. With the increase of time, the rates for all six segments of the wellbore reach a steady value.Figure 17
Flux distribution along the horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections (k1 = k3 = 10 mD).


## 3.3. Homogeneous-Dual-Porosity (Fracture-Matrix)-Triple-Porosity (Fracture-Vug-Matrix) Model

It is assumed that a horizontal well traverses a reservoir consisting of homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug- matrix) sections with the same permeability. The physical model is illustrated in Figure18.Figure 18
Homogeneous-dual-porosity (fracture-matrix)-triple-porosity (fracture-vug-matrix) model.In Figure19, the solid curves represent a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug-matrix) sections. The triangles represent the type curves for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections. When a horizontal well traverses a reservoir consisting of homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug- matrix) sections with equal permeability, due to the influence of the first homogeneous section and the middle dual-porosity section, the interporosity flow from vug into fracture becomes weaker than the interporosity flow from matrix into fracture. The dip of the interporosity flow in dual-porosity section becomes easier to recognize; however, the general flow regimes remain unchanged.Figure 19
Type curves of flow stage division for a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug-matrix) sections (solid curves). The triangles represent the type curves for a horizontal well traversing a three-section reservoir including homogeneous, triple-porosity (fracture-vug-matrix), and homogeneous sections.In the flux distribution in Figure20, due to the two interporosity flows of the triple-porosity section and the interporosity flow of the dual-porosity section in the early stage, the flow rate of the two segments of wellbore in the triple-porosity section and the flow rate of the two segments of wellbore in the dual-porosity section are larger than the flow rate of the two segments of wellbore in the homogeneous section. With the appearance of pseudoradial flow regime, the rates for all six segments of the wellbore become equal. At late time, the flux distribution tends to be stable and the segment closer to the heel of the horizontal well gets higher flux.Figure 20
Flow rate per segment for a horizontal well traverses homogeneous, dual-porosity (fracture-matrix), and triple-porosity (fracture-vug-matrix) sections.Based on this model, when the permeability of homogeneous section decreases from 100 mD to 40 mD and then to 10 mD, the values of pressure drop and pressure derivative increase (Figure21). In addition, the stage of the interporosity flow from both vug and matrix into fracture system becomes harder to recognize at early times. The dip of the interporosity flow in dual-porosity section becomes shallower. Other flow regimes remain unchanged.Figure 21
Type curves of flow stage division when the permeability of section 1 increases from 10 mD to 40 mD and then to 100 mD.When the permeability of dual-porosity section decreases from 200 mD to 100 mD and then to 50 mD, the values of pressure drop and pressure derivative increase as it can be observed in Figure22. In addition, the stage of the interporosity flow from both vug and matrix into fracture becomes harder to recognize at early times. However, the dip of the interporosity flow in dual-porosity section becomes easier to identify. Other flow regimes remain unchanged.Figure 22
Type curves of flow stage division when the permeability of section 2 increases from 50 mD to 100 mD and then to 200 mD.Figure23 shows that when the permeability of triple-porosity section decreases from 200 mD to 100 mD and then to 50 mD, the values of pressure drop and pressure derivative increase. In addition, the stage of the interporosity flow from both vug and matrix into fracture becomes easier to recognize ate early times. Moreover, the dips of the interporosity flows in triple-porosity section become deeper. Other flow regimes remain unchanged.Figure 23
Type curves of flow stage division when the permeability of section 3 increases from 50 mD to 100 mD and then to 200 mD.


## 4. Field Application

To show the advantage of the proposed model with a horizontal well traversing multisection formation with different properties, a field case in Tarim Basin is analyzed. The input parameters of the well are as follows. The length of the well is 590 m. Before the pressure buildup test, this well has produced about 150 days with an average daily production of 80 m3/d. Then, the pressure buildup test has been carried out for about 180 h.The log-log fitting curves are shown in Figure24. We select the model with a horizontal well traversing a three-section reservoir including homogeneous, dual-porosity, and homogeneous sections to fit the measured data. The permeability of the dual-porosity section is greater than that of homogeneous sections. We set the main parameters as x1 = 200 m, x2 = 160 m, and x3 = 300 m and the corresponding permeability of each section as K1 = 15 mD, K2 = 45 mD, and K3 = 15 mD to fit the data with the model. It can be seen that the agreement between the measured data and the fitting curve is excellent, except in the wellbore storage period. Because the manometers may be installed downward the middle of the formation, the wellbore storage period lasts for a short time and is not represented on the curve.Figure 24
The matched type curves of actual well test data.The main fitting results are as follows:KNf=45mD, Km = 15 mD, and λ = 2.66 × 10−7. From the fitting results and curves, it is shown that there is a dual-porosity section along the horizontal well. According to the microseismic monitoring results, it can be verified that there is a dual-porosity section in the formation consisting of fractures and pores.


## 5. Conclusions

A well test model for pressure transient analysis of horizontal wells extending through a carbonate reservoir consisting of natural fractures, rock matrix, and vugs with different properties was presented in this study. Considering the pseudosteady interporosity flows from rock matrix and vugs into fractures, a multizone triple-porosity model was established and solved by using the point source function, Green’s function, and coupling of multiple reservoir sections. The corresponding type curves were developed, and sensitivity analysis was carried out. This approach satisfies the need for pressure transient analysis for a horizontal well that traverses two or more regions with distinct properties in heterogeneous carbonate reservoirs. The following conclusions can be drawn:(1)
A horizontal well traversing a three-section reservoir including homogeneous and dual-porosity (fracture-matrix) sections identifies the stages of pseudosteady interporosity flow from matrix into fracture, fracture pseudoradial flow, system linear flow, system pseudoradial flow, and pseudosteady flow occur in sequence. The greater the difference of permeability between the dual-porosity section and the two homogeneous sections, the deeper the valley on the pressure derivative curve.(2)
A horizontal well traversing a three-section reservoir including homogeneous and triple-porosity (fracture-vug-matrix) sections recognizes the stages of pseudosteady interporosity flow from vug and matrix into fracture, fracture pseudoradial flow, system linear flow, system pseudoradial flow, and pseudosteady flow occur sequentially. The greater the permeability of the triple-porosity section, the deeper the dips of the interporosity flows in this section on the pressure derivative curve.(3)
It is shown that, in the early stage, the flow rate in the dual-porosity/triple-porosity section is larger than the flow rate in the homogeneous section. With the advent of the pseudoradial flow regime, the rates for all segments of the wellbore become equal. At late times, the flux distribution tends to be stable and the segment closer to the heel of the horizontal well gets higher flux.There are still some issues to be further explored. Acid fracturing techniques and hydraulic fracturing techniques have been widely practiced to enhance the well productivity and ultimate recovery for low permeability carbonate reservoirs and shale gas reservoirs. We can extend our model from the horizontal well to the multifractured horizontal well to study the effect of different fracture networks. Based on the pressure transient analysis of the multifractured horizontal well, we can further study the production decline analysis to improve the well productivity and ultimate recovery.



---

*Source: 3267458-2020-02-24.xml*