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

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 Figure 1 [33].

[figure omitted; refer to PDF]

$$\begin{array}{}\text{(1)}& \overline{G}\left(M_D,M_D^{\prime },s\right)=\frac{\mu }{2Kl{x}_{eD}{h}_D}\left\{\frac{ch\left(\sqrt{u}{\tilde{y}}_{D1}\right)+ch\left(\sqrt{u}{\tilde{y}}_{D2}\right)}{\sqrt{u}sh\left(\sqrt{u}{\tilde{y}}_{eD}\right)}+2{\displaystyle \sum _{k=1}^{\infty }\cos \left(\frac{k\pi x_D}{x_{eD}}\right)\cos \left(\frac{k\pi x_D^{\prime }}{x_{eD}}\right)\frac{ch\left({\epsilon }_k{\tilde{y}}_{D1}\right)+ch\left({\epsilon }_k{\tilde{y}}_{D2}\right)}{{\epsilon }_{k}sh\left(\sqrt{u}{\tilde{y}}_{eD}\right)}+2{\displaystyle \sum _{n=1}^{\infty }\cos \left(\frac{n\pi z_D}{z_{eD}}\right)\cos \left(\frac{n\pi z_D^{\prime }}{z_{eD}}\right)\left[\frac{ch\left({\epsilon }_n{\tilde{y}}_{D1}\right)+ch\left({\epsilon }_n{\tilde{y}}_{D2}\right)}{{\epsilon }_{n}sh\left(\sqrt{u}{\tilde{y}}_{eD}\right)}+2{\displaystyle \sum _{k=1}^{\infty }\cos \left(\frac{k\pi x_D}{x_{eD}}\right)\cos \left(\frac{k\pi x_D^{\prime }}{x_{eD}}\right)}\frac{ch\left({\epsilon }_{k,n}{\tilde{y}}_{D1}\right)+ch\left({\epsilon }_{k,n}{\tilde{y}}_{D2}\right)}{{\epsilon }_{k,n}sh\left(\sqrt{u}{\tilde{y}}_{eD}\right)}\right]}}\right\},\end{array}$$ Where$$\begin{array}{}\text{(2)}& h_{\text{D}}=z_{\text{eD}}=\frac{h}{l}\sqrt{\frac{K}{K_z}},\text{(3)}& {\tilde{y}}_{\text{D}1}=y_{\text{eD}}-\left|y_{\text{D}}-y_{\text{D}}^{\prime }\right|,\text{(4)}& {\tilde{y}}_{\text{D}2}=y_{\text{eD}}-\left|y_{\text{D}}+y_{\text{D}}^{\prime }\right|,\text{(5)}& {\epsilon }_n=\sqrt{u+{\left(\frac{n\pi }{h_{\text{D}}}\right)}^2},\text{(6)}& {\epsilon }_k=\sqrt{u+{\left(\frac{k\pi }{x_{\text{eD}}}\right)}^2},\text{(7)}& {\epsilon }_{k,n}=\sqrt{u+{\left(\frac{k\pi }{x_{\text{eD}}}\right)}^2+{\left(\frac{n\pi }{h_{\text{D}}}\right)}^2},\end{array}$$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, K_z 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$$\begin{array}{}\text{(8)}& u=\left\{\begin{array}{cc}s,\hfill & \text{for\hspace{0.17em}homogeneous\hspace{0.17em}reservoirs,}\hfill \\ sf\left(s\right),\hfill & \text{for\hspace{0.17em}naturally\hspace{0.17em}fractured\hspace{0.17em}reservoirs}\text{.}\hfill \end{array}\right.\end{array}$$

Appropriate expressions for $f\left(s\right)$ 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 $f\left(s\right)$ are given, respectively, as follows:$$\begin{array}{}\text{(9)}& f\left(s\right)=\frac{sw\left(1-w\right)+\lambda }{s\left(1-w\right)+\lambda },\text{(10)}& f\left(s\right)=w_{N\text{f}}+\frac{w_v{\lambda }_v}{s{w}_v+{\lambda }_v}+\frac{w_m{\lambda }_m}{s{w}_m+{\lambda }_m},\end{array}$$where $w$ is the dimensionless storativity and λ is the dimensionless transfer coefficient. The subscript N_f 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, B_w, is a line equal to the length of the horizontal well, L_h. Then, the source function in the Laplace domain, ${\overline{S}}_{\alpha }$, given by Medeiros and Ozkan can be expressed as$$\begin{array}{}\text{(11)}& {\overline{S}}_{wi}={\displaystyle {\int }_{x_{wi}-L_{hi}/2}^{x_{wi}+L_{hi}/2}\overline{G}\left(x_D,y_D,z_D,x_D^{\prime },y_{wiD},z_{wiD},s\right)\text{d}{x}^{\prime }},\end{array}$$where $x_{wi}$, $y_{wi}$, and $z_{wi}$ are the coordinates of the midpoint and $L_{hi}$ is the length of the i^th segment of the horizontal well. Substituting $\overline{G}$ given by equation (1) into equation (11) and evaluating the integral yield,$$\begin{array}{c}\text{(12)}& {\overline{S}}_{wi}\left(M_D,s\right)=C{L}_{hi}\frac{ch\left(\sqrt{u}{\tilde{y}}_{D1i}\right)+ch\left(\sqrt{u}{\tilde{y}}_{D2i}\right)}{\sqrt{u}sh\left(\sqrt{u}{\tilde{y}}_{eD}\right)}+2C{\displaystyle \sum _{k=1}^{\infty }\frac{ch\left({\epsilon }_k{\tilde{y}}_{D1i}\right)+ch\left({\epsilon }_k{\tilde{y}}_{D2i}\right)}{{\epsilon }_{k}sh\left({\epsilon }_k\sqrt{u}{\tilde{y}}_{eD}\right)}{\displaystyle {\int }_{x_{wi}-L_{hi}/2}^{x_{wi}+L_{hi}/2}\cos \left(\frac{k\pi x}{x_e}\right)\cos \left(\frac{k\pi x^{\prime }}{x_e}\right)\text{d}{x}^{\prime }}}\\ +2C{L}_{hi}{\displaystyle \sum _{n=1}^{\infty }\cos \left(\frac{n\pi z_D}{z_{eD}}\right)\cos \left(\frac{n\pi z_{wiD}}{z_{eD}}\right)\frac{ch\left({\epsilon }_n{\tilde{y}}_{D1i}\right)+ch\left({\epsilon }_n{\tilde{y}}_{D2i}\right)}{{\epsilon }_{n}sh\left({\epsilon }_n\sqrt{u}{\tilde{y}}_{eD}\right)}}\\ +4C{L}_{hi}{\displaystyle \sum _{n=1}^{\infty }\cos \left(\frac{n\pi z_D}{z_{eD}}\right)\cos \left(\frac{n\pi z_{wiD}}{z_{eD}}\right)}\\ \times {\displaystyle \sum _{j=1}^{\infty }\frac{ch\left({\epsilon }_{j,n}{\tilde{y}}_{D1i}\right)+ch\left({\epsilon }_{j,n}{\tilde{y}}_{D2i}\right)}{{\epsilon }_{j,n}sh\left({\epsilon }_{j,n}\sqrt{u}{\tilde{y}}_{eD}\right)}}{\displaystyle {\int }_{x_{wi}-L_{hi}/2}^{x_{wi}+L_{hi}/2}\cos \left(\frac{j\pi x}{x_e}\right)\cos \left(\frac{j\pi x^{\prime }}{x_e}\right)\text{d}{x}^{\prime }},\end{array}$$where$$\begin{array}{}\text{(13)}& C=\frac{\mu }{2Kl{x}_{eD}{h}_D}.\end{array}$$

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 the x-axis (a planar source in the y-z plane) is given by$$\begin{array}{}\text{(14)}& {\overline{S}}_{ej}={\displaystyle {\int }_{y_{ej}-y_{pj}/2}^{y_{ej}+y_{pj}/2}{\displaystyle {\int }_{z_{ej}-z_{pj}/2}^{z_{ej}+z_{pj}/2}\overline{G}\left(x_D,y_D,z_D,x_{wD},y^{\prime },z^{\prime },s\right)\text{d}{z}^{\prime }\text{d}{y}^{\prime }}},\end{array}$$where y_ej and z_ej are the coordinates of the midpoint and y_pj and z_pj are the lengths of the sides of the planar segment.

Substituting equation (1) into equation (14) and evaluating the internal result in$$\begin{array}{c}\text{(15)}& {\overline{S}}_{ej}\left(M_D,s\right)=C{z}_{pj}{\displaystyle {\int }_{y_{ej}-y_{pj}/2}^{y_{ej}+y_{pj}/2}\frac{ch\left(\sqrt{u}{\tilde{y}}_{D1j}\right)+ch\left(\sqrt{u}{\tilde{y}}_{D2j}\right)}{\sqrt{u}sh\left(\sqrt{u}{\tilde{y}}_{eD}\right)}\text{d}{y}^{\prime }}+2C{z}_{pj}{\displaystyle \sum _{k=1}^{\infty }\cos \left(\frac{k\pi x_D}{x_{eD}}\right)\cos \left(\frac{k\pi x_{wD}}{x_{eD}}\right){\displaystyle {\int }_{y_{ej}-y_{pj}/2}^{y_{ej}+y_{pj}/2}\frac{ch\left({\epsilon }_k{\tilde{y}}_{D1j}\right)+ch\left({\epsilon }_k{\tilde{y}}_{D2j}\right)}{{\epsilon }_{k}sh\left({\epsilon }_k{\tilde{y}}_{eD}\right)}\text{d}{y}^{\prime }}}\\ +2C{\displaystyle \sum _{n=1}^{\infty }{\displaystyle {\int }_{z_{ej}-z_{pj}/2}^{z_{ej}+z_{pj}/2}\cos \left(\frac{n\pi z}{z_e}\right)\cos \left(\frac{n\pi z^{\prime }}{z_e}\right)\text{d}{z}^{\prime }\times {\displaystyle {\int }_{y_{ej}-y_{pj}/2}^{y_{ej}+y_{pj}/2}\frac{ch\left({\epsilon }_n{\tilde{y}}_{D1j}\right)+ch\left({\epsilon }_n{\tilde{y}}_{D2j}\right)}{{\epsilon }_{n}sh\left({\epsilon }_n{\tilde{y}}_{eD}\right)}\text{d}{y}^{\prime }}}}\\ +4C{\displaystyle \underset{z_{ej}-z_{pj}/2}{\overset{z_{ej}+z_{pj}/2}{\int }}{\displaystyle \sum _{n=1}^{\infty }\cos \left(\frac{n\pi z}{z_e}\right)\cos \left(\frac{n\pi z^{\prime }}{z_e}\right)\text{d}{z}^{\prime }}}\\ \times {\displaystyle \sum _{k=1}^{\infty }\cos \left(\frac{k\pi x_D}{x_{eD}}\right)\cos \left(\frac{k\pi x_{wD}}{x_{eD}}\right){\displaystyle {\int }_{y_{ej}-y_{pj}/2}^{y_{ej}+y_{pj}/2}\frac{ch\left({\epsilon }_{k,n}{\tilde{y}}_{D1j}\right)+ch\left({\epsilon }_{k,n}{\tilde{y}}_{D2j}\right)}{{\epsilon }_{k,n}sh\left({\epsilon }_{k,n}{\tilde{y}}_{eD}\right)}\text{d}{y}^{\prime }}}.\end{array}$$

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 omitted; refer to PDF]

[figure omitted; refer to PDF]

From equation (1), the pressure drop at the center of a well or interface segment i in section k is given by$$\begin{array}{}\text{(16)}& \overline{\Delta p}\left(M_{Di},s\right)=\sum _{j=1}^2{\overline{\tilde{q}}}_{wj}^k\left(s\right){\overline{S}}_{iwj}^k\left(M_{Di},s\right)+{\displaystyle \sum _{l=1}^2}{\overline{\tilde{q}}}_{el}^k\left(s\right){\overline{S}}_{iel}^k\left(M{}_{Di},s\right),\end{array}$$where M_Di indicates the midpoint of segment i, ${\overline{\tilde{q}}}_{wj}$ and ${\overline{\tilde{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:$$\begin{array}{}\text{(17)}& {\overline{\tilde{q}}}_{w1}^k{\overline{S}}_{wiw1}^k+{\overline{\tilde{q}}}_{w2}^k{\overline{S}}_{wiw2}^k+{\overline{\tilde{q}}}_{e1}^k{\overline{S}}_{wie1}^k+{\overline{\tilde{q}}}_{e2}^k{\overline{S}}_{wie2}^k-{\overline{\Delta p}}_{wi}^k=0,\end{array}$$for i = 1, 2, k = 1, 3,$$\begin{array}{}\text{(18)}& {\overline{\tilde{q}}}_{w1}^k{\overline{S}}_{wiw1}^k+{\overline{\tilde{q}}}_{w2}^k{\overline{S}}_{wiw2}^k+{\overline{\tilde{q}}}_{e1}^{k1}{\overline{S}}_{wie1}^{k1}+{\overline{\tilde{q}}}_{e2}^{k1}{\overline{S}}_{wie2}^{k1}+{\overline{\tilde{q}}}_{e1}^{k2}{\overline{S}}_{wie1}^{k2}+{\overline{\tilde{q}}}_{e2}^{k2}{\overline{S}}_{wie2}^{k2}-{\overline{\Delta p}}_{wi}^k=0,\end{array}$$for i = 1, 2, k = 2,$$\begin{array}{}\text{(19)}& {\overline{\tilde{q}}}_{w1}^k{\overline{S}}_{eiw1}^k+{\overline{\tilde{q}}}_{w2}^k{\overline{S}}_{eiw2}^k+{\overline{\tilde{q}}}_{e1}^k{\overline{S}}_{eie1}^k+{\overline{\tilde{q}}}_{e2}^k{\overline{S}}_{eie2}^k-{\overline{\Delta p}}_{ei}^k=0,\end{array}$$for i = 1, 2, k = 1, 3,$$\begin{array}{}\text{(20)}& {\overline{\tilde{q}}}_{w1}^2{\overline{S}}_{eiw1}^{2k}+{\overline{\tilde{q}}}_{w2}^2{\overline{S}}_{eiw2}^{2k}+{\overline{\tilde{q}}}_{e1}^{21}{\overline{S}}_{eie1}^{2k1}+{\overline{\tilde{q}}}_{e2}^{21}{\overline{S}}_{eie2}^{2k1}-{\overline{\tilde{q}}}_{e1}^{22}{\overline{S}}_{eie1}^{2k2}+{\overline{\tilde{q}}}_{e2}^{22}{\overline{S}}_{eie2}^{2k2}-{\overline{\Delta p}}_{ei}^{2k}=0,\end{array}$$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$$\begin{array}{c}\text{(21)}& {\overline{\Delta p}}_{ei}^1={\overline{\Delta p}}_{ei}^{21},{\overline{\Delta p}}_{ei}^{22}={\overline{\Delta p}}_{ei}^3,\\ \text{for\hspace{0.17em}}i=\mathrm{1,2},\\ \text{(22)}& {\overline{\tilde{q}}}_{ei}^1=-{\overline{\tilde{q}}}_{ei}^{21},{\overline{\tilde{q}}}_{ei}^{22}=-{\overline{\tilde{q}}}_{ei}^3,\\ \text{for\hspace{0.17em}}i=\mathrm{1,2.}\end{array}$$

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:$$\begin{array}{c}\text{(23)}& {\overline{\Delta p}}_{w1}^k-{\overline{\Delta p}}_{w2}^k=0,\hspace{1em}\text{for\hspace{0.17em}}k=\mathrm{1,2,3},\text{(24)}& {\overline{\Delta p}}_{w2}^1-{\overline{\Delta p}}_{w1}^2=0,{\overline{\Delta p}}_{w2}^2-{\overline{\Delta p}}_{w1}^3=0.\end{array}$$

Finally, the mass balance requires that the sum of fluxes entering the wellbore is equal to the production rate at the sandface q:$$\begin{array}{}\text{(25)}& {\displaystyle \sum _{k=1}^3{\displaystyle \sum _{j=1}^2{\overline{\tilde{q}}}_{wj}^k{L}_{hj}=\frac{q}{s}}}.\end{array}$$

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 Figures 4 and 5, which includes n horizontal well segments and n − 1 interface segments between the sections, is used.

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

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:$$\begin{array}{}\text{(26)}& {\overline{\tilde{q}}}_{w1}^k{\overline{S}}_{wiw1}^k+{\overline{\tilde{q}}}_{w2}^k{\overline{S}}_{wiw2}^k+{\overline{\tilde{q}}}_{e1}^k{\overline{S}}_{wie1}^k+{\overline{\tilde{q}}}_{e2}^k{\overline{S}}_{wie2}^k-{\overline{\Delta p}}_{wi}^k=0,\end{array}$$for i = 1, 2, k = 1, n,$$\begin{array}{}\text{(27)}& {\overline{\tilde{q}}}_{w1}^k{\overline{S}}_{wiw1}^k+{\overline{\tilde{q}}}_{w2}^k{\overline{S}}_{wiw2}^k+{\overline{\tilde{q}}}_{e1}^{k1}{\overline{S}}_{wie1}^{k1}+{\overline{\tilde{q}}}_{e2}^{k1}{\overline{S}}_{wie2}^{k1}+{\overline{\tilde{q}}}_{e1}^{k2}{\overline{S}}_{wie1}^{k2}+{\overline{\tilde{q}}}_{e2}^{k2}{\overline{S}}_{wie2}^{k2}-{\overline{\Delta p}}_{wi}^k=0,\end{array}$$for i = 1, 2, k = 2, 3, …, n − 1,$$\begin{array}{}\text{(28)}& {\overline{\tilde{q}}}_{w1}^k{\overline{S}}_{eiw1}^k+{\overline{\tilde{q}}}_{w2}^k{\overline{S}}_{eiw2}^k+{\overline{\tilde{q}}}_{e1}^k{\overline{S}}_{eie1}^k+{\overline{\tilde{q}}}_{e2}^k{\overline{S}}_{eie2}^k-{\overline{\Delta p}}_{ei}^k=0,\end{array}$$for i = 1, 2, k = 1, n,$$\begin{array}{}\text{(29)}& {\overline{\tilde{q}}}_{w1}^k{\overline{S}}_{eiw1}^{km}+{\overline{\tilde{q}}}_{w2}^k{\overline{S}}_{eiw2}^{km}+{\overline{\tilde{q}}}_{e1}^{k1}{\overline{S}}_{eie1}^{km1}+{\overline{\tilde{q}}}_{e2}^{k1}{\overline{S}}_{eie2}^{km1}-{\overline{\tilde{q}}}_{e1}^{k2}{\overline{S}}_{eie1}^{km2}+{\overline{\tilde{q}}}_{e2}^{k2}{\overline{S}}_{eie2}^{km2}-{\overline{\Delta p}}_{ei}^{km}=0,\end{array}$$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$$\begin{array}{c}\text{(30)}& {\overline{\Delta p}}_{ei}^1={\overline{\Delta p}}_{ei}^{21},{\overline{\Delta p}}_{ei}^{k2}={\overline{\Delta p}}_{ei}^{\left(k+1\right)1},\\ {\overline{\Delta p}}_{ei}^{\left(n-1\right)2}={\overline{\Delta p}}_{ei}^n,\\ \text{for\hspace{0.17em}}i=\mathrm{1,2},k=\mathrm{2,3},\dots ,n-1,\\ \text{(31)}& {\overline{\tilde{q}}}_{ei}^1=-{\overline{\tilde{q}}}_{ei}^{21},{\overline{\tilde{q}}}_{ei}^{k2}=-{\overline{\tilde{q}}}_{ei}^{\left(k+1\right)1},\\ {\overline{\tilde{q}}}_{ei}^{\left(n-1\right)2}=-{\overline{\tilde{q}}}_{ei}^n,\\ \text{for\hspace{0.17em}}i=\mathrm{1,2},k=\mathrm{2,3},\dots ,n-1.\end{array}$$

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:$$\begin{array}{}\text{(32)}& {\overline{\Delta p}}_{w1}^k-{\overline{\Delta p}}_{w2}^k=0,\hspace{1em}\text{for\hspace{0.17em}}k=\mathrm{1,2,3},\dots ,n,\text{(33)}& {\overline{\Delta p}}_{w2}^k-{\overline{\Delta p}}_{w1}^{k+1}=0,\hspace{1em}\text{for\hspace{0.17em}}k=\mathrm{1,2,3},\dots ,n-1.\end{array}$$

Finally, mass balance requires that the sum of fluxes entering the wellbore is equal to the production rate at the sandface q:$$\begin{array}{}\text{(34)}& {\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^2{\overline{\tilde{q}}}_{wj}^k{L}_{hj}=\frac{q}{s}}}.\end{array}$$

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 = {A₁, A₂}, are given by$$\begin{array}{c}\text{(35)}& A_2=\left[\begin{array}{cccccc}-1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & -1\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & -1\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -1\hfill \\ 1\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 1\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill & -1\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 1\hfill & -1\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill & -1\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right],A_1=\left[\begin{array}{cccccccccc}{\overline{S}}_{w1w1}^1\hfill & {\overline{S}}_{w1w2}^1\hfill & {\overline{S}}_{w1e1}^1\hfill & {\overline{S}}_{w1e2}^1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ {\overline{S}}_{w2w1}^1\hfill & {\overline{S}}_{w2w2}^1\hfill & {\overline{S}}_{w2e1}^1\hfill & {\overline{S}}_{w2e2}^1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ {\overline{S}}_{e1w1}^1\hfill & {\overline{S}}_{e1w2}^1\hfill & -\left({\overline{S}}_{e1e1}^1+{\overline{S}}_{e1e1}^{211}\right)\hfill & -\left({\overline{S}}_{e1e2}^1+{\overline{S}}_{e1e2}^{211}\right)\hfill & -{\overline{S}}_{e1w1}^{21}\hfill & -{\overline{S}}_{e1w2}^{21}\hfill & {\overline{S}}_{e1e1}^{212}\hfill & {\overline{S}}_{e1e2}^{212}\hfill & 0\hfill & 0\hfill \\ {\overline{S}}_{e2w1}^1\hfill & {\overline{S}}_{e2w2}^1\hfill & -\left({\overline{S}}_{e2e1}^1+{\overline{S}}_{e2e1}^{211}\right)\hfill & -\left({\overline{S}}_{e2e2}^1+{\overline{S}}_{e2e2}^{211}\right)\hfill & -{\overline{S}}_{e2w1}^{21}\hfill & -{\overline{S}}_{e2w2}^{21}\hfill & {\overline{S}}_{e2e1}^{212}\hfill & {\overline{S}}_{e2e2}^{212}\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & {\overline{S}}_{w1e1}^{21}\hfill & {\overline{S}}_{w1e2}^{21}\hfill & {\overline{S}}_{w1w1}^2\hfill & {\overline{S}}_{w1w2}^2\hfill & -{\overline{S}}_{w1e1}^{22}\hfill & -{\overline{S}}_{w1e2}^{22}\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & {\overline{S}}_{w2e1}^{21}\hfill & {\overline{S}}_{w2e2}^{21}\hfill & {\overline{S}}_{w2w1}^2\hfill & {\overline{S}}_{w2w2}^2\hfill & -{\overline{S}}_{w2e1}^{22}\hfill & -{\overline{S}}_{w2e2}^{22}\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & {\overline{S}}_{e1e1}^{221}\hfill & {\overline{S}}_{e1e2}^{221}\hfill & {\overline{S}}_{e1w1}^{22}\hfill & {\overline{S}}_{e1w2}^{22}\hfill & -\left({\overline{S}}_{e1e1}^{222}+{\overline{S}}_{e1e1}^3\right)\hfill & -\left({\overline{S}}_{e1e2}^{222}+{\overline{S}}_{e1e2}^3\right)\hfill & {\overline{S}}_{e1w1}^3\hfill & {\overline{S}}_{e1w2}^3\hfill \\ 0\hfill & 0\hfill & {\overline{S}}_{e2e1}^{221}\hfill & {\overline{S}}_{e2e2}^{221}\hfill & {\overline{S}}_{e2w1}^{22}\hfill & {\overline{S}}_{e1w2}^{22}\hfill & -\left({\overline{S}}_{e2e1}^{222}+{\overline{S}}_{e2e1}^3\right)\hfill & -\left({\overline{S}}_{e2e2}^{222}+{\overline{S}}_{e2e2}^3\right)\hfill & {\overline{S}}_{e2w1}^3\hfill & {\overline{S}}_{e2w2}^3\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & {\overline{S}}_{w1e1}^3\hfill & {\overline{S}}_{w1e2}^3\hfill & {\overline{S}}_{w1w1}^3\hfill & {\overline{S}}_{w1w2}^3\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & {\overline{S}}_{w2e1}^3\hfill & {\overline{S}}_{w2e2}^3\hfill & {\overline{S}}_{w2w1}^3\hfill & {\overline{S}}_{w2w2}^3\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ L_{h1}^1\hfill & L_{h2}^1\hfill & 0\hfill & 0\hfill & L_{h1}^2\hfill & L_{h2}^2\hfill & 0\hfill & 0\hfill & L_{h1}^3\hfill & L_{h2}^3\hfill \end{array}\right].\end{array}$$

The solution vector, x = [x₁, x₂], has the following components:$$\begin{array}{c}\text{(36)}& x_1=\left[\begin{array}{cccccccccc}{\overline{\tilde{q}}}_{w1}^1\hfill & {\overline{\tilde{q}}}_{w2}^1\hfill & {\overline{\tilde{q}}}_{e1}^{21}\hfill & {\overline{\tilde{q}}}_{e2}^{21}\hfill & {\overline{\tilde{q}}}_{w1}^2\hfill & {\overline{\tilde{q}}}_{w1}^2\hfill & {\overline{\tilde{q}}}_{e1}^3\hfill & {\overline{\tilde{q}}}_{e2}^3\hfill & {\overline{\tilde{q}}}_{w1}^3\hfill & {\overline{\tilde{q}}}_{w1}^3\hfill \end{array}\right],x_2=\left[\begin{array}{cccccc}{\overline{\Delta p}}_{w1}^1\hfill & {\overline{\Delta p}}_{w2}^1\hfill & {\overline{\Delta p}}_{w1}^2\hfill & {\overline{\Delta p}}_{w2}^2\hfill & {\overline{\Delta p}}_{w1}^3\hfill & {\overline{\Delta p}}_{w2}^3\hfill \end{array}\right].\end{array}$$

The components of the right-hand-side vector, b = [b₁, b₂], are$$\begin{array}{c}\text{(37)}& b_1=\left[\begin{array}{cccccccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right],b_2=\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & q/s\hfill \end{array}\right].\end{array}$$

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 = {A₁, A₂}, are given by$$\begin{array}{c}\text{(38)}& A_2=\left[\begin{array}{cccccccccc}-1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & -1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & 0\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & 0\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & -1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & -1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & 0\hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & 0\hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \ddots \hfill & 4\left(n-3\right)\text{rows}\hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & -1\hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & -1\hfill & \hfill \\ 1\hfill & -1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & 1\hfill & -1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & 1\hfill & -1\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \ddots \hfill & \ddots \hfill & 2\left(n-3\right)\text{rows}\hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 1\hfill & -1\hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 1\hfill & -1\hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 1\hfill & -1\hfill \\ 0\hfill & \cdots \hfill & \hfill & \left(2n-2\text{\hspace{0.17em}columns}\right)\hfill & \hfill & \hfill & \hfill & \hfill & \cdots \hfill & 0\hfill \end{array}\right],A_1=\left[\begin{array}{cccccccccccc}\hfill & {\overline{S}}_{w1w1}^{\left(k-1\right)}\hfill & {\overline{S}}_{w1w2}^{\left(k-1\right)}\hfill & {\overline{S}}_{w1e1}^{\left(k-1\right)2}\hfill & {\overline{S}}_{w1e2}^{\left(k-1\right)2}\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & {\overline{S}}_{w2w1}^{\left(k-1\right)}\hfill & {\overline{S}}_{w2w2}^{\left(k-1\right)}\hfill & {\overline{S}}_{w2e1}^{\left(k-1\right)2}\hfill & {\overline{S}}_{w2e2}^{\left(k-1\right)2}\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & {\overline{S}}_{e1w1}^{\left(k-1\right)}\hfill & {\overline{S}}_{e1w2}^{\left(k-1\right)2}\hfill & -\left({\overline{S}}_{e1e1}^{\left(k-1\right)22}+{\overline{S}}_{e1e1}^{k11}\right)\hfill & -\left({\overline{S}}_{e1e2}^{\left(k-1\right)22}+{\overline{S}}_{e1e2}^{k11}\right)\hfill & -{\overline{S}}_{e1w1}^{k1}\hfill & -{\overline{S}}_{e1w2}^{k1}\hfill & {\overline{S}}_{e1e1}^{k12}\hfill & {\overline{S}}_{e1e2}^{k12}\hfill & \hfill & \hfill & \hfill \\ \hfill & {\overline{S}}_{e2w1}^{\left(k-1\right)}\hfill & {\overline{S}}_{e2w2}^{\left(k-1\right)2}\hfill & -\left({\overline{S}}_{e2e1}^{\left(k-1\right)22}+{\overline{S}}_{e2e1}^{k11}\right)\hfill & -\left({\overline{S}}_{e2e2}^{\left(k-1\right)22}+{\overline{S}}_{e2e2}^{k11}\right)\hfill & -{\overline{S}}_{e2w1}^{k1}\hfill & -{\overline{S}}_{e2w2}^{k1}\hfill & {\overline{S}}_{e2e1}^{k12}\hfill & {\overline{S}}_{e2e2}^{k12}\hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & {\overline{S}}_{w1e1}^{k1}\hfill & {\overline{S}}_{w1e2}^{k1}\hfill & {\overline{S}}_{w1w1}^k\hfill & {\overline{S}}_{w1w2}^k\hfill & -{\overline{S}}_{w1e1}^{k2}\hfill & -{\overline{S}}_{w1e2}^{k2}\hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & {\overline{S}}_{w2e1}^{k1}\hfill & {\overline{S}}_{w2e2}^{k1}\hfill & {\overline{S}}_{w2w1}^k\hfill & {\overline{S}}_{w2w2}^k\hfill & -{\overline{S}}_{w2e1}^{k2}\hfill & -{\overline{S}}_{w2e2}^{k2}\hfill & \hfill & \hfill & \hfill \\ \cdots \hfill & \hfill & \hfill & {\overline{S}}_{e1e1}^{k21}\hfill & {\overline{S}}_{e1e2}^{k21}\hfill & {\overline{S}}_{e1w1}^{k2}\hfill & {\overline{S}}_{e1w2}^{k2}\hfill & -\left({\overline{S}}_{e1e1}^{k22}+{\overline{S}}_{e1e1}^{\left(k+1\right)11}\right)\hfill & -\left({\overline{S}}_{e1e2}^{k22}+{\overline{S}}_{e1e2}^{\left(k+1\right)11}\right)\hfill & {\overline{S}}_{e1w1}^{\left(k+1\right)1}\hfill & {\overline{S}}_{e1w2}^{\left(k+1\right)1}\hfill & \cdots \hfill \\ \hfill & \hfill & \hfill & {\overline{S}}_{e2e1}^{k21}\hfill & {\overline{S}}_{e2e2}^{k21}\hfill & {\overline{S}}_{e2w1}^{k2}\hfill & {\overline{S}}_{e2w2}^{k2}\hfill & -\left({\overline{S}}_{e2e1}^{k22}+{\overline{S}}_{e2e1}^{\left(k+1\right)11}\right)\hfill & -\left({\overline{S}}_{e2e2}^{k22}+{\overline{S}}_{e2e2}^{\left(k+1\right)11}\right)\hfill & {\overline{S}}_{e2w1}^{\left(k+1\right)1}\hfill & {\overline{S}}_{e2w2}^{\left(k+1\right)1}\hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & {\overline{S}}_{w1e1}^{\left(k+1\right)1}\hfill & {\overline{S}}_{w1e2}^{\left(k+1\right)1}\hfill & {\overline{S}}_{w1w1}^{\left(k+1\right)}\hfill & {\overline{S}}_{w1w2}^{\left(k+1\right)}\hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & {\overline{S}}_{w2e1}^{\left(k+1\right)1}\hfill & {\overline{S}}_{w2e2}^{\left(k+1\right)1}\hfill & {\overline{S}}_{w2w1}^{\left(k+1\right)}\hfill & {\overline{S}}_{w2w2}^{\left(k+1\right)}\hfill & \hfill \\ \hfill & 0\hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \ddots \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill \\ \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & \hfill & 0\hfill & \hfill \\ \cdots \hfill & L_{h1}^{k-1}\hfill & L_{h2}^{k-1}\hfill & 0\hfill & 0\hfill & L_{h1}^k\hfill & L_{h2}^k\hfill & 0\hfill & 0\hfill & L_{h1}^{k+1}\hfill & L_{h2}^{k+1}\hfill & \cdots \hfill \end{array}\right].\end{array}$$

The solution vector, x = [x₁, x₂], has the following components:$$\begin{array}{c}\text{(39)}& x_1=\left[\begin{array}{cccccccccccc}{\overline{\tilde{q}}}_{w1}^1\hfill & {\overline{\tilde{q}}}_{w2}^1\hfill & \cdots \hfill & {\overline{\tilde{q}}}_{e1}^{k1}\hfill & {\overline{\tilde{q}}}_{e2}^{k1}\hfill & {\overline{\tilde{q}}}_{w1}^k\hfill & {\overline{\tilde{q}}}_{w1}^k\hfill & \cdots \hfill & {\overline{\tilde{q}}}_{e1}^n\hfill & {\overline{\tilde{q}}}_{e2}^n\hfill & {\overline{\tilde{q}}}_{w1}^n\hfill & {\overline{\tilde{q}}}_{w1}^n\hfill \end{array}\right],\hspace{1em}\text{for\hspace{0.17em}}k=\mathrm{2,3},\dots ,n-1,x_2=\left[\begin{array}{cccccccc}{\overline{\Delta p}}_{w1}^1\hfill & {\overline{\Delta p}}_{w2}^1\hfill & \cdots \hfill & {\overline{\Delta p}}_{w1}^k\hfill & {\overline{\Delta p}}_{w2}^k\hfill & \cdots \hfill & {\overline{\Delta p}}_{w1}^n\hfill & {\overline{\Delta p}}_{w2}^n\hfill \end{array}\right],\hspace{1em}\text{for\hspace{0.17em}}k=\mathrm{2,3},\dots ,n-1.\end{array}$$

The components of the right-hand-side vector, b = [b₁, b₂], are$$\begin{array}{c}\text{(40)}& b_1={\left[\begin{array}{cccccccc}0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill \end{array}\right]}_{1\times 4n-2},b_2={\left[\begin{array}{cccccc}0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & q/s\hfill \end{array}\right]}_{1\times 2n}.\end{array}$$

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 Figure 6.

[figure omitted; refer to PDF]

In Figure 7, 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 omitted; refer to PDF]

It can be seen from Figure 8 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 omitted; refer to PDF]

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 Figure 9.

[figure omitted; refer to PDF]

In Figure 9, 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 Figure 10, 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 omitted; refer to PDF]

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 Figure 11. 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 omitted; refer to PDF]

Figure 12 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 omitted; refer to PDF]

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 Figure 13.

[figure omitted; refer to PDF]

In Figure 14, 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 omitted; refer to PDF]

It can be observed from the flux distribution in Figure 15 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 omitted; refer to PDF]

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 Figure 16.

[figure omitted; refer to PDF]

In Figure 16, 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 Figure 17 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 omitted; refer to PDF]

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 Figure 18.

[figure omitted; refer to PDF]

In Figure 19, 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 omitted; refer to PDF]

In the flux distribution in Figure 20, 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 omitted; refer to PDF]

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 (Figure 21). 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 omitted; refer to PDF]

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 Figure 22. 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 omitted; refer to PDF]

Figure 23 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 omitted; refer to PDF]

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 m³/d. Then, the pressure buildup test has been carried out for about 180 h.

The log-log fitting curves are shown in Figure 24. 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 x₁ = 200 m, x₂ = 160 m, and x₃ = 300 m and the corresponding permeability of each section as K₁ = 15 mD, K₂ = 45 mD, and K₃ = 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 omitted; refer to PDF]

The main fitting results are as follows: $K_{N_{\text{f}}}=45\hspace{0.17em}\text{mD}$, K_m = 15 mD, and λ = 2.66 × 10⁻⁷. 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.

Acknowledgments

This work was supported by a grant from the National Science and Technology Major Project of the Ministry of Science and Technology of China (no. 2016ZX05015-003-004-05).

Glossary

Nomenclature

μ:Fluid viscosity (mPa·s)

K:Permeability (10⁻³ μm²)

K_z:Permeability in the z direction (10⁻³μm²)

l:Reference length (m)

L_hi:Length of the i^th segment (m)

Δp:Pressure drop (MPa)

q:Production rate at the sandface (m³/d)

$q_w$:Production of wellbore (m³/d)

q _e:Production of outer boundary (m³/d)

M_Di:Midpoint of segment i

$w$:Fluid storage capacitance coefficient

λ:Dimensionless transfer coefficient

x:x coordinate (m)

y:y coordinate (m)

z:z coordinate (m)

s:Laplace transform parameter

Δ:Difference operator

N_f:Natural fracture system

$v$:Vug system

m:Matrix system

D:Dimensionless

e:Boundary

$w$:Wellbore

n:Number

j:Number

k:Number

l:Number.

References

[1] T. Gentzis, D. Bolen, "The use of numerical simulation in predicting coalbed methane producibility from the Gates coals, Alberta Inner Foothills, Canada: comparison with Mannville CBM production in Alberta Syncline," International Journal of Coal Geology, vol. 74 no. 3-4, pp. 215-236, DOI: 10.1016/j.coal.2007.12.003, 2008.

[2] T. Gentzis, N. Deisman, R. J. Chalaturnyk, "A method to predict geomechanical properties and model well stability in horizontal boreholes," International Journal of Coal Geology, vol. 78 no. 2, pp. 149-160, DOI: 10.1016/j.coal.2008.11.001, 2009.

[3] I. Palmer, "Coalbed methane: a world view," International Journal of Coal Geology, vol. 82 no. 3-4, pp. 184-195, DOI: 10.1016/j.coal.2009.12.010, 2010.

[4] Z. Pan, L. D. Connell, "Modelling permeability for coal reservoirs: a review of analytical models and testing data," International Journal of Coal Geology, vol. 92,DOI: 10.1016/j.coal.2011.12.009, 2012.

[5] M. E. Houben, N. J. Hardebol, A. Barnhoorn, "Fluid flow from matrix to fractures in Early Jurassic shales," International Journal of Coal Geology, vol. 175, pp. 26-39, DOI: 10.1016/j.coal.2017.03.012, 2017.

[6] G. I. Barenblatt, I. P. Kochina, P. Zheltov, "Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]," Journal of Applied Mathematics and Mechanics, vol. 24 no. 5, pp. 1286-1303, DOI: 10.1016/0021-8928(60)90107-6, 1960.

[7] J. E. Warren, P. J. Root, "The behavior of naturally fractured reservoirs," Society of Petroleum Engineers Journal, vol. 3 no. 3, pp. 245-255, DOI: 10.2118/426-pa, 1963.

[8] H. Kazemi, "Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution," Society of Petroleum Engineers Journal, vol. 9 no. 4, pp. 451-462, DOI: 10.2118/2156-a, 1969.

[9] A. O. de Swaan, "Analytic solutions for determining naturally fractured reservoir properties by well testing," Society of Petroleum Engineers Journal, vol. 16 no. 3, pp. 117-122, DOI: 10.2118/5346-pa, 1976.

[10] C. Braester, "Influence of block size on the transition curve for a drawdown test in a naturally fractured reservoir," Society of Petroleum Engineers Journal, vol. 24 no. 5, pp. 498-504, DOI: 10.2118/10543-pa, 1984.

[11] H. Cinco-Ley, F. Samaniego, F. Kucuk, "The pressure transient behavior for naturally fractured reservoirs with multiple block size," Proceedings of the SPE Annual Technical Conference and Exhibition,DOI: 10.2118/14168-ms, .

[12] A. K. Belani, Y. Jalali, "Estimation of matrix block size distribution in naturally fractured reservoirs," Proceedings of the SPE Annual Technical Conference and Exhibition,DOI: 10.2118/18171-ms, .

[13] J. Spivey, W. Lee, "Pressure transient response for a naturally fractured reservoir with a distribution of block sizes," Proceedings of the SPE Rocky Mountain Regional/Low-Permeability Reservoirs Symposium and Exhibition,DOI: 10.2118/60294-ms, .

[14] R. Raghavan, C. U. Ohaeri, "Unsteady flow to a well produced at a constant pressure in a fractured reservoir," SPE Formation Evaluation, vol. 2 no. 2, pp. 186-200, 1981.

[15] R. C. Velazquez, M. A. Vasquez, R. Castrejonr-Aivar, V. Arana-Ortiz, "Pressure- transient and decline-curve behavior in naturally fractured vuggy carbonate reservoirs," SPE Reservoir Evaluation and Engineering, vol. 8 no. 2, pp. 95-111, DOI: 10.2118/77689-pa, 2005.

[16] D. Abdassah, I. Ershaghi, "Triple-porosity systems for representing naturally fractured reservoirs," SPE Formation Evaluation, vol. 1 no. 2, pp. 113-127, DOI: 10.2118/13409-pa, 1986.

[17] Y. Jalali, I. Ershaghi, "Pressure transient analysis of heterogeneous naturally fractured reservoirs," Proceedings of the SPE California Regional Meeting,DOI: 10.2118/16341-ms, .

[18] A. Al-Ghamdi, I. Ershaghi, "Pressure transient analysis of dually fractured reservoirs," SPE Journal, vol. 1 no. 1, pp. 93-100, DOI: 10.2118/26959-pa, 1996.

[19] C. R. Clarkson, "Production data analysis of unconventional gas wells: Workflow," International Journal of Coal Geology, vol. 109-110, pp. 147-157, DOI: 10.1016/j.coal.2012.11.016, 2013.

[20] J. Ren, L. Zhang, S. Ren, "Multi-branched horizontal wells for coalbed methane production: field performance and well structure analysis," International Journal of Coal Geology, vol. 131, pp. 52-64, DOI: 10.1016/j.coal.2014.06.003, 2014.

[21] F. Medeiros, E. Ozkan, H. Kazemi, "A semianalytical, pressure-transient model for horizontal and multilateral wells in composite, layered, and compartmentalized reservoirs," Proceedings of the SPE Annual Technical Conference and Exhibition,DOI: 10.2118/102834-ms, .

[22] C. H. Bowman, P. B. Crawford, A Practical Method for Calculating the Transient Pressure Distribution in Linear Semi-infinite Water Driven Reservoirs Having Different Fluid and Rock Properties in Each Zone, 1961.

[23] H. C. Bixel, B. K. Larkin, H. K. Van Poollen, "Effect of linear discontinuities on pressure build-up and drawdown behavior," Journal of Petroleum Technology, vol. 15 no. 8, pp. 885-895, DOI: 10.2118/611-pa, 1963.

[24] T. D. Strltsova, R. M. Mckinley, "Effect of flow time duration on buildup pattern for reservoirs with heterogeneous properties," Society of Petroleum Engineers Journal, vol. 24 no. 3, pp. 294-306, DOI: 10.2118/11140-pa, 1984.

[25] A. K. Ambastha, P. G. Mcleroy, A. S. Grader, "Effects of a partially communicating fault in a composite reservoir on transient pressure testing," SPE Formation Evaluation, vol. 4 no. 2, pp. 210-218, DOI: 10.2118/16764-pa, 1989.

[26] M. J. Bourgeois, F. H. Daviau, J.-L. Boutaud de la Combe, "Pressure behavior in finite channel-levee complexes," SPE Formation Evaluation, vol. 11 no. 3, pp. 177-184, DOI: 10.2118/26461-pa, 1996.

[27] F. J. Kuchuk, H. Tarek, "Pressure behavior of laterally composite reservoirs," SPE Formation Evaluation, vol. 12 no. 12, pp. 47-56, DOI: 10.2118/24678-pa, 1997.

[28] M. Wei, Y. Duan, M. Dong, Q. Fang, "Blasingame decline type curves with material balance pseudo-time modified for multi-fractured horizontal wells in shale gas reservoirs," Journal of Natural Gas Science and Engineering, vol. 31, pp. 340-350, DOI: 10.1016/j.jngse.2016.03.033, 2016.

[29] M. Wei, Y. Duan, Q. Fang, T. Zhang, "Production decline analysis for a multi-fractured horizontal well considering elliptical reservoir stimulated volumes in shale gas reservoirs," Journal of Geophysics and Engineering, vol. 13 no. 3, pp. 354-365, DOI: 10.1088/1742-2132/13/3/354, 2016.

[30] M. Dejam, H. Hassanzadeh, Z. Chen, "Semi-analytical solutions for a partially penetrated well with wellbore storage and skin effects in a double-porosity system with a gas cap," Transport in Porous Media, vol. 100 no. 2, pp. 159-192, DOI: 10.1007/s11242-013-0210-6, 2013.

[31] M. Dejam, H. Hassanzadeh, Z. Chen, "Semi-analytical solution for pressure transient analysis of a hydraulically fractured vertical well in a bounded dual-porosity reservoir," Journal of Hydrology, vol. 565, pp. 289-301, DOI: 10.1016/j.jhydrol.2018.08.020, 2018.

[32] M. Dejam, H. Hassanzadeh, "The role of natural fractures of finite double-porosity aquifers on diffusive leakage of brine during geological storage of CO 2," International Journal of Greenhouse Gas Control, vol. 78, pp. 177-197, DOI: 10.1016/j.ijggc.2018.08.007, 2018.

[33] E. Ozkan, R. Raghavan, "New solutions for well-test-analysis problems: Part 1-analytical considerations (includes associated papers 28666 and 29213)," SPE Formation Evaluation, vol. 6 no. 3, pp. 359-368, DOI: 10.2118/18615-pa, 1991.

[34] E. Ozkan, R. Raghavan, "New solutions for well-test-analysis problems: part 2 computational considerations and applications," SPE Formation Evaluation, vol. 6 no. 3, pp. 369-378, DOI: 10.2118/18616-pa, 1991.

[35] R. Raghavan, Well Test Analysis Prentice Hall Petroleum Engineering Series, 1993.

[36] H. Stehfest, "Algorithm 368: numerical inversion of Laplace transforms [D5]," Communications of the Acm, vol. 13 no. 1, pp. 47-49, DOI: 10.1145/361953.361969, 1970.