# Analysis of Gas Flow Behavior for Highly Deviated Wells in Naturally Fractured-Vuggy Carbonate Gas Reservoirs

Authors: Kongjie Wang; Lian Wang; Caspar Daniel Adenutsi; Zhiping Li; Sen Yang; Liang Zhang; Lan Wang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

To improve the carbonate gas reservoir development and production, highly deviated wells (HDW) are widely used in the field. Production decline analysis of HDW is crucial for long-term gas reservoir development. However, it is a new challenge to incorporate the complex pore structure of naturally fractured-vuggy carbonate gas reservoirs and evaluate the production performance of HDW. This paper presents a semianalytical model to analyze the pressure and production behavior of HDW in naturally fractured-vuggy carbonate gas reservoirs, which consist of fractures, vugs, and matrix. The primary flow occurs only through the fracture and the outer boundary is closed. Introducing pseudopressure and pseudotime, the Laplace transformation, Fourier transformation, and its inverse and Stehfest numerical inversion were employed to establish a point source and line source solutions. Furthermore, the validity of the proposed model was verified by comparing a field data from the Arum River Basin in Turkmenistan. Finally, the effects of major parameters on the production decline curves were analyzed by using the proposed model and it was found that they had influences at different stages of gas production history and the sensitivity intensity of each parameter was different. With its high efficiency and simplicity, this semianalytical model will serve as a useful tool to evaluate the well production behavior for the naturally fractured-vuggy carbonate gas reservoirs.


---

## Body

## 1. Introduction

Carbonate gas reservoirs are widely distributed throughout the world. In these carbonate reservoirs, the original gas in place (OGIP) and gas production account for more than half of gas reserves and production in the world [1]. Carbonate gas reservoirs commonly have extremely complex porous structures. This mainly comprises three types of pores, which are matrix, natural fractures, and vugs with different degrees of development. Thus, gas flow in carbonate gas reservoirs show characteristics of multiple-porosity system [2]. To improve gas reservoir development and production, highly deviated wells (HDW) are widely used [3]. A new challenge that arises from the application of this technology to improve reservoir development is the difficulty in analyzing and predicting the production behaviors of HDW in naturally fractured-vuggy carbonate gas reservoirs.In the past decades, a lot of researchers reported some relevant works. The theory and applications on dual-porosity flow model for fractured reservoirs have been well researched [4–8]. The dual-porosity model was first proposed by Warren and Root in which it was assumed that the reservoir has matrix blocks with low permeability but high porosity and the fracture system with high permeability but low porosity [9]. Based on this model, Bourdet and Gringarten introduced a log-log analysis model with wellbore storage and skin in dual-porosity systems [10]. Jalali and Ershaghi established a pressure transient model for carbonate reservoirs on semilog plots and found that the slope characteristics of the transition cover an embracing spectrum [11]. Ei-Banbi extended previous analytical models to develop a catalog of linear flow causes, models, and solutions for dual-porosity linear reservoirs, dual-porosity radial reservoirs, and dual-porosity channel reservoirs [12]. Ying-Lan et al. represented a dual-porosity model within a porous-vuggy carbonate reservoir that does not introduce the fracture system [13].Reservoirs composed of matrix, fracture, and vugs are called triple-porosity systems. A lot of triple-porosity models have been proposed through analytical, semianalytical, and numerical methods for analyzing gas flow in fractured-vuggy carbonate reservoirs. With regard to analytical methods, a triple-porosity and single-permeability model was first proposed by Abdassah and Ershaghi [14]. The authors considered an unsteady-state interporosity flow model between the fractures and the others, which is a more realistic representation. Al-Ghamdi and Ershaghi presented a triple-porosity dual-permeability model to characterize fracture heterogeneities, which considered discontinuous matrix and two continuous fracture networks, microfractures and macrofractures [15]. Later, another triple-porosity dual-permeability model was proposed to characterize vuggy porosity, coexisting with matrix and fractures in carbonate reservoir by Camacho-Velázquez et al. [16]. Wu et al. presented an analytical method for pressure transient analysis in fractured-vuggy carbonate reservoir that relied on a triple-continuum concept [17]. Lately, Wang et al. developed a model for transient flow analysis of acid fracturing wells in fractured-vuggy carbonate reservoirs [18]. With regard to semianalytical and numerical methods, Gulbransen et al. proposed a multiscale mixed finite element method for modeling fractured-vuggy carbonate reservoirs [19]. Li et al. built a coupled Stoked-Darcy model for estimating the equivalent permeability tensor and modeling two-phase flow with a full permeability tensor in a fracture-vug- media system [20]. Zhang et al. presented a new numerical model that considered the effect of geomechanics on fractures and vugs [21]. For upscaling fractured-vuggy reservoir models, Gao et al. introduced a numerical model that consisted of four different porosity systems, i.e., the matrix, fractures, isolated vugs, and connected vugs [22].A substantial amount of research has focused on the pressure behavior and production performance of deviated wells. Abbaszadeh and Hegeman first proposed analytical solution for the pressure drawdown behavior of a deviated limited-entry well with different combinations of boundaries. The authors used the method of source and Green’s functions to derive these solutions [23]. Ozkan and Raghavan presented a solution for a deviated limited-entry well in an infinite reservoir with no-flow top and bottom boundaries by employing the Laplace transformation [24]. Harmohan et al. introduced a method of numerical simulation to analyze the pressure transient behavior of deviated wells that intersect high permeability layers between two low permeability layers [25]. Based on rigorous derivation, Wang et al. presented a transient pressure solution for deviated wells, which reduced the amount of computation and significantly improved the computational efficiency [26]. In a recent work, Meng et al. proposed a semianalytical model of HDW for evaluating the pressure behavior which considered the stress-sensitive performance in fractured-vuggy carbonate gas reservoirs with composite systems [27]. Zhao et al. built a new model for evaluating inflow performance of HDW in anisotropic reservoirs, which was proposed by introducing the anisotropic quantitative characterization method [28].Currently, only few suitable analytical models exist to evaluate the production behavior of HDW in fractured-vuggy carbonate gas reservoirs. In this study, a semianalytical flow model was established by employing Laplace transformation, Fourier transformation and inversion, Stehfest numerical inversion algorithm, and point source function. Based on the proposed model, the effects of major parameters on production performance were studied.


## 2. Physical Model

Figure1 shows a schematic for HDW in a fractured-vuggy carbonate gas reservoir. Fractured-vuggy carbonate gas reservoirs are naturally structured by matrix system, natural fractures system, and vugs system. Figure 2 shows the physical modeling scheme of a fractured-vuggy carbonate medium. In order to mathematically define the gas flow behavior, the assumptions of this model are listed as follows:(1)
HDW are located in a laterally infinite gas reservoir.(2)
The gas reservoir is assumed to be horizontal with equal thickness and have two impermeable boundaries at the top and bottom.(3)
During the well production, gas flows only through the fracture system into the wellbore; interporosity flow occurs from the vugs and the matrix systems to the fracture system.(4)
Flow is single phase, isothermal and follows Darcy's law. Gravitational and capillary effects are negligible.(5)
The pressure is uniform initially throughout the gas reservoir and is equal topi.(6)
The HDW produces at a constant rateqg or at constant pressure pwf.(7)
The permeabilitykh in horizontal direction and kz in vertical direction are different.Figure 1
Schematic diagram of HDW in a fractured-vuggy carbonate gas reservoir.Figure 2
Flow sketch map for a cell [29].


## 3. Mathematical Model

### 3.1. Establishment of Mathematical Model

In this work, the interporosity flow from vugs and matrix to fractures is assumed to be pseudosteady state. Combining the continuity equation with transport equation as well as the equation of state and the introduction of pseudopressure and pseudotime (the detailed procedure for pseudopressure and pseudotime is documented in AppendixA), the governing equation for the triple-porosity system in radial cylindrical coordinate is given as follows:(1)1r∂∂rr∂mf∂r+kvfkf∂2mf∂z2=10.0864ϕfμgCtfkhf∂mf∂tp-αmkmkhmm-mf-αckckhmc-mfFor the matrix system the equation is(2)ϕmμgCtm0.0864∂mm∂tp=αmkmmf-mmFor the vugs system the equation is(3)ϕcμgCtc0.0864∂mc∂tp=αckcmf-mcwheremf, mm, mc is the pseudopressure of fractures system, matrix system, vugs system, respectively, MPa/s; ϕf, ϕm, ϕc is porosity of fracture system, matrix system, vugs system, respectively, fraction; kf, km, kc is permeability of fracture system, matrix system, vugs system, respectively,md; kvf is vertical fracture permeability,md; Ctf, Ctm, Ctc is compressibility of fracture system, matrix system, vugs system, respectively, MPa-1; μg is gas viscosity, mPa·s; αm, αc are shape factors of matrix and vugs, 1/m2; and tp is pseudotime, day.The initial pseudopressuremi is assumed uniform for the triple-porosity system:(4)mjr,t=0=mi,j=f,m,cThe outer boundary is closed:(5)∂mj∂rr=re=0,j=f,m,cAt the inner boundary, there is a continuous point source (xw,yw,zw) with a constant gas rate q~:(6)q~Bgi=limε→0limr→0khfαpμgε∫zw-ε/2zw+ε/2r∂mf∂rdzz-zw≤ε20z-zw>ε2wherekhf is fracture horizontal permeability, md.At the top and bottom, the boundaries are impermeable:(7)∂mf∂zz=0=∂mf∂zz=h=0whereh is formation thickness, m.Based on the definitions of dimensionless variables in Table1, the governing equations and the initial and boundary conditions (Eq. (1)-(7)) can be transformed into dimensionless terms (Eq. (8)-(14)).Table 1
Definitions of dimensionless variables.
Dimensionless variables Definitions Dimensionless pseudopressure m j D = π k f h T s c m i ( p i ) - m j ( p ) p s c q g T  ( j = f , m , c )  Dimensionless pseudotime t D = α t k f t p μ g r w 2 ( ϕ C t ) f + m + c  Dimensionless radius r D = r r w k f k h f  Dimensionless formation radius r e D = r e r w k f k h f  Dimensionless distance ofx coordinate x D = x r w k f k h f  Dimensionless distance ofy coordinate y D = y r w k f k h f  Dimensionless distance of mid-perforation inx coordinate x w D = x w r w k f k h f  Dimensionless distance of mid-perforation iny coordinate y w D = y w r w k f k h f  Dimensionless formation thickness h D = h r w k f k h f  Dimensionless vertical distance z D = z r w k f k h f  Dimensionless vertical distance of mid-perforation z w D = z w r w k f k h f  Dimensionless infinitesimal variable ε D = ε r w k f k h f  Dimensionless deviated well length L w D = L w r w k f k h f s i n 2 θ + k f k v f c o s 2 θ  Dimensionless storativity ratio ω j = ϕ j C t j ( ϕ C t ) f + m + c ,  ( j = f , m , c )  Dimensionless interporosity flow coefficient λ j = α j k j r w 2 k f ,  ( j = m , c ) Dimensionless production rate q w D = p s c T π T s c k f h m i ( p i ) - m w ( p w f ) q gThe dimensionless governing equation of fracture, matrix, and vugs systems are(8)1rD∂∂rDrD∂mfD∂rD+∂2mfD∂zD2=ωf∂mfD∂tD-λmmmD-mfD-λcmcD-mfD(9)ωm∂mmD∂tD=λmmfD-mmD(10)ωc∂mcD∂tD=λcmfD-mcDTransformed initial condition:(11)mjDrD,tD=0=0,j=f,m,cTransformed outer boundary condition:(12)∂mjD∂rDrD=reD=0,j=f,m,cTransformed inner boundary condition:(13)limεD→0⁡limrD→0⁡∫zwD-εD/2zwD+εD/2rD∂mfD∂rDdzD=-hDzD-zwD≤εD20zD-zwD>εD2Transformed top and bottom boundaries:(14)∂mfD∂zDzD=0=∂mfD∂zDzD=hD=0


### 3.2. Solution of Mathematical Model

#### 3.2.1. Basic Continuous Point Source Solution

Laplace transformation, Fourier transformation, and inversion were employed to solve the dimensionless model (Eq. (8)-(14)). The solution in Laplace space is as follows (the detailed procedure for the solution is documented in Appendix B).(15)m¯fD=1sK1reDsfsI1reDsfsI0rDsfs+K0rDsfs+2s∑n=1∞K1reDun2+sfsI1reDun2+sfsI0rDun2+sfs+K0rDun2+sfscos⁡μnzwDcos⁡μnzDwhere(16)un=nπhD(17)fs=wf+λcωcλc+ωcs+λmωmλm+ωmswherek0 is modified Bessel function of zero order of the second kind; k1 is modified Bessel function of order one of the second kind; I0 is modified Bessel function of order zero of the first kind; I1 is modified Bessel function of order one of the first kind.


#### 3.2.2. Pressure Distribution of a HDW

To obtain the pressure solution of HDWs, the wellbore was treated as a uniform flux line source. It was assumed that there was an infinitesimal point on the HDW. Based on the principle of superposition for a point source, integration was carried out along the deviated line for the line source solution. Cinco (1975) successfully got ther, θ, z-coordinates of the equivalent-pressure point where the wellbore pressure of HDW can be calculated easily [30]. Therefore, the dimensionless pressure distribution of HDW in Laplace domain is given as follows:(18)m¯fD=1sLwD∫-LwD/2LwD/2K1reDsfsI1reDsfsI0r~Dsfs+K0r~Dsfsdη+2sLwD∫-LwD/2LwD/2∑n=1∞K1reDun2+sfsI1reDun2+sfsI0r~Dun2+sfs+K0r~Dun2+sfscosμnz~wDcosμnzDdηwhere(19)r~D=xD-xwD-ηsin⁡θw2+yD-ywD2(20)z~wD=zwD+ηcos⁡θw(21)θw=tan-1⁡kvfkhftan⁡θ(22)LwD=Lwrwkfkhfsin2θ+kfkvfcos2θBased on the principle of superposition, the pressure solution with wellbore storage and skin effect in Laplace domain can be obtained as(23)m¯wD=sm¯fD+Ss+CDs2sm¯fD+S(24)CD=C6.28ϕCtf+c+mhrw2whereC is wellbore storage coefficient, m3/MPa; m¯wD is the dimensionless bottom pseudopressure.


#### 3.2.3. Gas Production Based on the Well Pressure Responses

According to Van Everdingen and Hurst [31], the dimensionless production response for a HDW producing at a constant bottom-hole pressure in the Laplace domain can be obtained as follows:(25)q¯wD=1s2m¯wDEquation (25) can then be inverted to obtain qwD in real space using suitable numerical inversion algorithms such as Stehfest's inversion algorithm [32].


## 3.1. Establishment of Mathematical Model

In this work, the interporosity flow from vugs and matrix to fractures is assumed to be pseudosteady state. Combining the continuity equation with transport equation as well as the equation of state and the introduction of pseudopressure and pseudotime (the detailed procedure for pseudopressure and pseudotime is documented in AppendixA), the governing equation for the triple-porosity system in radial cylindrical coordinate is given as follows:(1)1r∂∂rr∂mf∂r+kvfkf∂2mf∂z2=10.0864ϕfμgCtfkhf∂mf∂tp-αmkmkhmm-mf-αckckhmc-mfFor the matrix system the equation is(2)ϕmμgCtm0.0864∂mm∂tp=αmkmmf-mmFor the vugs system the equation is(3)ϕcμgCtc0.0864∂mc∂tp=αckcmf-mcwheremf, mm, mc is the pseudopressure of fractures system, matrix system, vugs system, respectively, MPa/s; ϕf, ϕm, ϕc is porosity of fracture system, matrix system, vugs system, respectively, fraction; kf, km, kc is permeability of fracture system, matrix system, vugs system, respectively,md; kvf is vertical fracture permeability,md; Ctf, Ctm, Ctc is compressibility of fracture system, matrix system, vugs system, respectively, MPa-1; μg is gas viscosity, mPa·s; αm, αc are shape factors of matrix and vugs, 1/m2; and tp is pseudotime, day.The initial pseudopressuremi is assumed uniform for the triple-porosity system:(4)mjr,t=0=mi,j=f,m,cThe outer boundary is closed:(5)∂mj∂rr=re=0,j=f,m,cAt the inner boundary, there is a continuous point source (xw,yw,zw) with a constant gas rate q~:(6)q~Bgi=limε→0limr→0khfαpμgε∫zw-ε/2zw+ε/2r∂mf∂rdzz-zw≤ε20z-zw>ε2wherekhf is fracture horizontal permeability, md.At the top and bottom, the boundaries are impermeable:(7)∂mf∂zz=0=∂mf∂zz=h=0whereh is formation thickness, m.Based on the definitions of dimensionless variables in Table1, the governing equations and the initial and boundary conditions (Eq. (1)-(7)) can be transformed into dimensionless terms (Eq. (8)-(14)).Table 1
Definitions of dimensionless variables.
Dimensionless variables Definitions Dimensionless pseudopressure m j D = π k f h T s c m i ( p i ) - m j ( p ) p s c q g T  ( j = f , m , c )  Dimensionless pseudotime t D = α t k f t p μ g r w 2 ( ϕ C t ) f + m + c  Dimensionless radius r D = r r w k f k h f  Dimensionless formation radius r e D = r e r w k f k h f  Dimensionless distance ofx coordinate x D = x r w k f k h f  Dimensionless distance ofy coordinate y D = y r w k f k h f  Dimensionless distance of mid-perforation inx coordinate x w D = x w r w k f k h f  Dimensionless distance of mid-perforation iny coordinate y w D = y w r w k f k h f  Dimensionless formation thickness h D = h r w k f k h f  Dimensionless vertical distance z D = z r w k f k h f  Dimensionless vertical distance of mid-perforation z w D = z w r w k f k h f  Dimensionless infinitesimal variable ε D = ε r w k f k h f  Dimensionless deviated well length L w D = L w r w k f k h f s i n 2 θ + k f k v f c o s 2 θ  Dimensionless storativity ratio ω j = ϕ j C t j ( ϕ C t ) f + m + c ,  ( j = f , m , c )  Dimensionless interporosity flow coefficient λ j = α j k j r w 2 k f ,  ( j = m , c ) Dimensionless production rate q w D = p s c T π T s c k f h m i ( p i ) - m w ( p w f ) q gThe dimensionless governing equation of fracture, matrix, and vugs systems are(8)1rD∂∂rDrD∂mfD∂rD+∂2mfD∂zD2=ωf∂mfD∂tD-λmmmD-mfD-λcmcD-mfD(9)ωm∂mmD∂tD=λmmfD-mmD(10)ωc∂mcD∂tD=λcmfD-mcDTransformed initial condition:(11)mjDrD,tD=0=0,j=f,m,cTransformed outer boundary condition:(12)∂mjD∂rDrD=reD=0,j=f,m,cTransformed inner boundary condition:(13)limεD→0⁡limrD→0⁡∫zwD-εD/2zwD+εD/2rD∂mfD∂rDdzD=-hDzD-zwD≤εD20zD-zwD>εD2Transformed top and bottom boundaries:(14)∂mfD∂zDzD=0=∂mfD∂zDzD=hD=0


## 3.2. Solution of Mathematical Model

### 3.2.1. Basic Continuous Point Source Solution

Laplace transformation, Fourier transformation, and inversion were employed to solve the dimensionless model (Eq. (8)-(14)). The solution in Laplace space is as follows (the detailed procedure for the solution is documented in Appendix B).(15)m¯fD=1sK1reDsfsI1reDsfsI0rDsfs+K0rDsfs+2s∑n=1∞K1reDun2+sfsI1reDun2+sfsI0rDun2+sfs+K0rDun2+sfscos⁡μnzwDcos⁡μnzDwhere(16)un=nπhD(17)fs=wf+λcωcλc+ωcs+λmωmλm+ωmswherek0 is modified Bessel function of zero order of the second kind; k1 is modified Bessel function of order one of the second kind; I0 is modified Bessel function of order zero of the first kind; I1 is modified Bessel function of order one of the first kind.


### 3.2.2. Pressure Distribution of a HDW

To obtain the pressure solution of HDWs, the wellbore was treated as a uniform flux line source. It was assumed that there was an infinitesimal point on the HDW. Based on the principle of superposition for a point source, integration was carried out along the deviated line for the line source solution. Cinco (1975) successfully got ther, θ, z-coordinates of the equivalent-pressure point where the wellbore pressure of HDW can be calculated easily [30]. Therefore, the dimensionless pressure distribution of HDW in Laplace domain is given as follows:(18)m¯fD=1sLwD∫-LwD/2LwD/2K1reDsfsI1reDsfsI0r~Dsfs+K0r~Dsfsdη+2sLwD∫-LwD/2LwD/2∑n=1∞K1reDun2+sfsI1reDun2+sfsI0r~Dun2+sfs+K0r~Dun2+sfscosμnz~wDcosμnzDdηwhere(19)r~D=xD-xwD-ηsin⁡θw2+yD-ywD2(20)z~wD=zwD+ηcos⁡θw(21)θw=tan-1⁡kvfkhftan⁡θ(22)LwD=Lwrwkfkhfsin2θ+kfkvfcos2θBased on the principle of superposition, the pressure solution with wellbore storage and skin effect in Laplace domain can be obtained as(23)m¯wD=sm¯fD+Ss+CDs2sm¯fD+S(24)CD=C6.28ϕCtf+c+mhrw2whereC is wellbore storage coefficient, m3/MPa; m¯wD is the dimensionless bottom pseudopressure.


### 3.2.3. Gas Production Based on the Well Pressure Responses

According to Van Everdingen and Hurst [31], the dimensionless production response for a HDW producing at a constant bottom-hole pressure in the Laplace domain can be obtained as follows:(25)q¯wD=1s2m¯wDEquation (25) can then be inverted to obtain qwD in real space using suitable numerical inversion algorithms such as Stehfest's inversion algorithm [32].


## 3.2.1. Basic Continuous Point Source Solution

Laplace transformation, Fourier transformation, and inversion were employed to solve the dimensionless model (Eq. (8)-(14)). The solution in Laplace space is as follows (the detailed procedure for the solution is documented in Appendix B).(15)m¯fD=1sK1reDsfsI1reDsfsI0rDsfs+K0rDsfs+2s∑n=1∞K1reDun2+sfsI1reDun2+sfsI0rDun2+sfs+K0rDun2+sfscos⁡μnzwDcos⁡μnzDwhere(16)un=nπhD(17)fs=wf+λcωcλc+ωcs+λmωmλm+ωmswherek0 is modified Bessel function of zero order of the second kind; k1 is modified Bessel function of order one of the second kind; I0 is modified Bessel function of order zero of the first kind; I1 is modified Bessel function of order one of the first kind.


## 3.2.2. Pressure Distribution of a HDW

To obtain the pressure solution of HDWs, the wellbore was treated as a uniform flux line source. It was assumed that there was an infinitesimal point on the HDW. Based on the principle of superposition for a point source, integration was carried out along the deviated line for the line source solution. Cinco (1975) successfully got ther, θ, z-coordinates of the equivalent-pressure point where the wellbore pressure of HDW can be calculated easily [30]. Therefore, the dimensionless pressure distribution of HDW in Laplace domain is given as follows:(18)m¯fD=1sLwD∫-LwD/2LwD/2K1reDsfsI1reDsfsI0r~Dsfs+K0r~Dsfsdη+2sLwD∫-LwD/2LwD/2∑n=1∞K1reDun2+sfsI1reDun2+sfsI0r~Dun2+sfs+K0r~Dun2+sfscosμnz~wDcosμnzDdηwhere(19)r~D=xD-xwD-ηsin⁡θw2+yD-ywD2(20)z~wD=zwD+ηcos⁡θw(21)θw=tan-1⁡kvfkhftan⁡θ(22)LwD=Lwrwkfkhfsin2θ+kfkvfcos2θBased on the principle of superposition, the pressure solution with wellbore storage and skin effect in Laplace domain can be obtained as(23)m¯wD=sm¯fD+Ss+CDs2sm¯fD+S(24)CD=C6.28ϕCtf+c+mhrw2whereC is wellbore storage coefficient, m3/MPa; m¯wD is the dimensionless bottom pseudopressure.


## 3.2.3. Gas Production Based on the Well Pressure Responses

According to Van Everdingen and Hurst [31], the dimensionless production response for a HDW producing at a constant bottom-hole pressure in the Laplace domain can be obtained as follows:(25)q¯wD=1s2m¯wDEquation (25) can then be inverted to obtain qwD in real space using suitable numerical inversion algorithms such as Stehfest's inversion algorithm [32].


## 4. Results and Discussion

### 4.1. Proposed Model Application and Validation

Field data from the Arum River Basin in Turkmenistan was used to validate the proposed model and to demonstrate its practical application. Arum River Basin is a large scale sedimentary basin, which is located between the border of Turkmenistan and Uzbekistan. The reservoir bed possesses carbonate rock that is formed by postdepositional diagenesis, including dissolution and dolomitization. Vugs, fractures, and dissolution pores are highly dispersed in the reservoir. The details of relevant parameters needed for the model evaluation are listed in Table2. These parameters are obtained by various methods in the field. In particular, it should be noted that formation radius (re), fracture permeability (kf), interporosity flow coefficients between different media systems (λc,λm), and initial reservoir pressure (pi) can be obtained through interpretation of well testing data. And the porosity of three different media (ϕf,ϕm,ϕc) can be obtained by combining well logging data with core CT scanning.Table 2
The values of relevant parameters for verification and sensitivity analysis.
Parameters Value Wellbore radius,rw (m) 0.0889 Formation radius,re (m) 400 Formation thickness,h(m) 30 Fracture permeability,kf (m) 1 Interporosity flow coefficient between fractures and vugs,λc 2 × 1 0 - 4 Interporosity flow coefficient between fractures and matrix,λm 2 × 1 0 - 6 Fractures porosity,ϕf 0.002 Matrix porosity,ϕm 0.04 Vugs porosity,ϕc 0.017 Inclination angle,θ (°) 75 Length of the highly deviated well,Lw (m) 300 Initial reservoir pressure,pi (MPa) 24 Bottom-hole pressure,pw (MPa) 20 Reservoir temperature,T (K) 381.15 Gas production rate,qg (m3/d) 600,000 Gas viscosity,μg (mPa·s) 0.02 Gas compressibility factor,Z (m3/m3) 0.97Figure3 shows a comparison between daily rate of gas production from a highly deviated well from the field and simulated gas production rate from the proposed model. It was found that the proposed model had a good fitting with actual production data. Except for very few data points that deviated from the simulated curve due to the real production time being less than 24 hours, the rest of the data fits well with the curve. After ignoring data points less than 24 hours of production, the average relative error (E) between the remaining actual data and model data is only 4.19%. The average relative error (E) is evaluated by(26)E=1n∑j=1nqgj-qgj~qgj~Figure 3
Comparing simulated results and production rates of Arum River Basin HDW.whereqgi is simulated production rate from the proposed model, m3/d; qgj~ is field production rate, m3/d.Moreover, it could be seen that the trend of simulated production rate in the following 1000 days and even longer predicted the gas production trend for the future with a constant bottom-hole pressure. Based on this application, it is confirmed that the proposed model can be used to describe the gas production behavior of HDW in naturally fractured-vuggy carbonate gas reservoir. Therefore, the estimated input parameters of the model will be used as the basis for the sensitivity analysis.


### 4.2. Modeling of Production Behavior with Triple-Porosity Flows

In order to understand the production behavior with triple-porosity flows, different flow regime analysis was conducted to study the effects of various flow mechanisms on the overall production behavior. Figure4 shows the complete transient flow process of HDW production in fractured-vuggy carbonate gas reservoir under closed circle boundary and it can be divided into five flow stages. The first stage was dominated by wellbore storage and skin effect and the dimensionless pseudopressure and pseudopressure derivative curves showed a straight line with a slope of 1. After that, skin factor affected the shape of the derivative curve that looked like a “hump”. The second stage was dominated by inclination angle (θ) of HDW. When the inclination angle approached 0°, the HDW was treated as a vertical well. In this case, this stage would not be obvious or would completely disappear on the curve. The third stage was dominated by interporosity flow between fractures and vugs and this flow regime behaved with a “dip” on the pressure derivative curve. Because the gas flow in the vugs was relatively smoother than in the matrix, the interporosity flow between fractures and vugs appeared first when the wellbore pressure began to deplete. Therefore, wellbore pressure depletion was retarded due to gas supplement from vugs to fracture system. The fourth stage was dominated by interporosity flow between fractures and matrix and this stage also behaved with a “dip” on the pressure derivative curve. It is worth noting that the third and fourth stages could interfere with each other, depending on the value of λc and λm. The fifth stage was dominated by closed circle boundary, which manifested on the derivative curve with a slope of 1.Figure 4
Typical pressure curves for HDW in fractured-vuggy carbonate gas reservoir under closed circle boundary.


### 4.3. Model Sensitivity Analyses

Based on the field data, the production sensitivity analyses for different key parameters are discussed. In this section, the HDW produces with a constant bottom-hole pressure. The key parameters include formation thickness, fracture permeability, inclination angle of the HDW, formation radius, interporosity flow coefficient between fractures and vugs and vugs storativity ratio. The influences of these parameters are apparent and their estimates are essential for future production analysis and forecasting.Figure5 shows the effects of formation thickness on production performance. From Figure 5, it can be seen that formation thickness primarily influences every stage of gas production. For the time, t<100 days, the difference of gas production rates between different formation thicknesses were almost the same. It meant that the decline rates of gas production for different scenarios were almost the same. After the time, t>100 days, the differences become smaller with increase in time. In addition, it was also observed that with decrease in formation thickness, the area enclosed by the curve became smaller, which meant the cumulative production reduced.Figure 5
Effect of formation thickness on production decline curves.From Figure6, it was noted that permeability affects the whole period of gas production. The higher the fracture permeability, the higher the initial production rate. However, for a reservoir with higher fracture permeability, production rate would decrease drastically, and during the middle-late stages the production rate would be exceeded by a reservoir with lower fracture permeability. Moreover, a well with low fracture permeability has a long stable production period.Figure 6
Effect of fracture permeability on production decline curves.Figure7 shows the effects of reservoir radius on gas production performance of HDWs. The difference of gas production at late stage was increasingly wider. It was noted that with increase in gas reservoir radius, the production curves declined more gently, which meant a gas reservoir with bigger radius would have a longer stable production period.Figure 7
Effect of gas reservoir radius on production decline curves.Figure8 shows the effects of interporosity flow coefficient between fractures and vugs on gas production performance of the HDW. It can be seen that the coefficient primarily affected early stage of gas production. For the time, t<10 days, with the interporosity flow coefficient between fractures and vugs smaller, the initial production also became smaller and the production rate decreased relatively slower. For the time, t>10 days, there was no difference between the production rates of the three scenarios. From Figure 9, it was observed that vugs storativity ratio affected middle-late stage of gas production. During 10<t<10000 days, with vugs storativity ratio smaller, production rate decreased more rapidly. However, after t>10000 days, the difference in the production curves was contrary to the previous one. Meanwhile, it was obvious that the extent of production variation is less drastic than the middle stage.Figure 8
Effect of interporosity flow coefficient between fractures and vugs on production decline curves.Figure 9
Effect of vugs storativity ratio on production decline curves.Figure10 shows the effects of the angle of inclination on gas production performance of the HDW. It can be seen that the angle of inclination also affects the entire stage of gas production. Due to the effects of different inclination angle on contact area of wellbore, the cumulative gas production increases with increase in inclination angle.Figure 10
Effect of inclination angle on production decline curves.According to Figures5–10, the production decline curve was affected by many parameters and the intensity of production variation was different at different stages of gas production history. By extracting production rate data at different times (t=10days, 100days, 1000days) for adjacent parameters in Table 3, production difference percentage (Δqper) could be obtained which could be helpful in comparing the sensitivity of these parameters.Table 3
Effect of parameters att=10 days, 100 days, 1000 days on the difference in production rate percentage.
parameters  10 days 100 days 1000 days Δ q p e r, % Δ q p e r, % Δ q p e r, % k f, md 1/0.5 22.1 18.4 18.2 r e, m 800/400 3.8 12.8 20.9 λ c 0.0002/0.00002 4.6 1.1 1.2 ω c 0.1/0.01 3.7 7.5 5.9 θ, degree 80/40 23.1 21.8 18.7Based on Table3, trends of difference in percentage of parameters are shown in Figure 11. It can be seen that the bigger the difference in percentage value, the more sensitive the parameters. Inclination angle (θ) and fracture permeability (kf) seemed to be the two most sensitive parameters; reservoir radius (re) was the third, vugs storativity ratio (ωc) and interporosity flow coefficient between fractures and vugs (λc) were the least two sensitive parameters.Figure 11
Trends of percentage difference of parameters by extracting production rate att=10 days, 100 days, 1000 days.


## 4.1. Proposed Model Application and Validation

Field data from the Arum River Basin in Turkmenistan was used to validate the proposed model and to demonstrate its practical application. Arum River Basin is a large scale sedimentary basin, which is located between the border of Turkmenistan and Uzbekistan. The reservoir bed possesses carbonate rock that is formed by postdepositional diagenesis, including dissolution and dolomitization. Vugs, fractures, and dissolution pores are highly dispersed in the reservoir. The details of relevant parameters needed for the model evaluation are listed in Table2. These parameters are obtained by various methods in the field. In particular, it should be noted that formation radius (re), fracture permeability (kf), interporosity flow coefficients between different media systems (λc,λm), and initial reservoir pressure (pi) can be obtained through interpretation of well testing data. And the porosity of three different media (ϕf,ϕm,ϕc) can be obtained by combining well logging data with core CT scanning.Table 2
The values of relevant parameters for verification and sensitivity analysis.
Parameters Value Wellbore radius,rw (m) 0.0889 Formation radius,re (m) 400 Formation thickness,h(m) 30 Fracture permeability,kf (m) 1 Interporosity flow coefficient between fractures and vugs,λc 2 × 1 0 - 4 Interporosity flow coefficient between fractures and matrix,λm 2 × 1 0 - 6 Fractures porosity,ϕf 0.002 Matrix porosity,ϕm 0.04 Vugs porosity,ϕc 0.017 Inclination angle,θ (°) 75 Length of the highly deviated well,Lw (m) 300 Initial reservoir pressure,pi (MPa) 24 Bottom-hole pressure,pw (MPa) 20 Reservoir temperature,T (K) 381.15 Gas production rate,qg (m3/d) 600,000 Gas viscosity,μg (mPa·s) 0.02 Gas compressibility factor,Z (m3/m3) 0.97Figure3 shows a comparison between daily rate of gas production from a highly deviated well from the field and simulated gas production rate from the proposed model. It was found that the proposed model had a good fitting with actual production data. Except for very few data points that deviated from the simulated curve due to the real production time being less than 24 hours, the rest of the data fits well with the curve. After ignoring data points less than 24 hours of production, the average relative error (E) between the remaining actual data and model data is only 4.19%. The average relative error (E) is evaluated by(26)E=1n∑j=1nqgj-qgj~qgj~Figure 3
Comparing simulated results and production rates of Arum River Basin HDW.whereqgi is simulated production rate from the proposed model, m3/d; qgj~ is field production rate, m3/d.Moreover, it could be seen that the trend of simulated production rate in the following 1000 days and even longer predicted the gas production trend for the future with a constant bottom-hole pressure. Based on this application, it is confirmed that the proposed model can be used to describe the gas production behavior of HDW in naturally fractured-vuggy carbonate gas reservoir. Therefore, the estimated input parameters of the model will be used as the basis for the sensitivity analysis.


## 4.2. Modeling of Production Behavior with Triple-Porosity Flows

In order to understand the production behavior with triple-porosity flows, different flow regime analysis was conducted to study the effects of various flow mechanisms on the overall production behavior. Figure4 shows the complete transient flow process of HDW production in fractured-vuggy carbonate gas reservoir under closed circle boundary and it can be divided into five flow stages. The first stage was dominated by wellbore storage and skin effect and the dimensionless pseudopressure and pseudopressure derivative curves showed a straight line with a slope of 1. After that, skin factor affected the shape of the derivative curve that looked like a “hump”. The second stage was dominated by inclination angle (θ) of HDW. When the inclination angle approached 0°, the HDW was treated as a vertical well. In this case, this stage would not be obvious or would completely disappear on the curve. The third stage was dominated by interporosity flow between fractures and vugs and this flow regime behaved with a “dip” on the pressure derivative curve. Because the gas flow in the vugs was relatively smoother than in the matrix, the interporosity flow between fractures and vugs appeared first when the wellbore pressure began to deplete. Therefore, wellbore pressure depletion was retarded due to gas supplement from vugs to fracture system. The fourth stage was dominated by interporosity flow between fractures and matrix and this stage also behaved with a “dip” on the pressure derivative curve. It is worth noting that the third and fourth stages could interfere with each other, depending on the value of λc and λm. The fifth stage was dominated by closed circle boundary, which manifested on the derivative curve with a slope of 1.Figure 4
Typical pressure curves for HDW in fractured-vuggy carbonate gas reservoir under closed circle boundary.


## 4.3. Model Sensitivity Analyses

Based on the field data, the production sensitivity analyses for different key parameters are discussed. In this section, the HDW produces with a constant bottom-hole pressure. The key parameters include formation thickness, fracture permeability, inclination angle of the HDW, formation radius, interporosity flow coefficient between fractures and vugs and vugs storativity ratio. The influences of these parameters are apparent and their estimates are essential for future production analysis and forecasting.Figure5 shows the effects of formation thickness on production performance. From Figure 5, it can be seen that formation thickness primarily influences every stage of gas production. For the time, t<100 days, the difference of gas production rates between different formation thicknesses were almost the same. It meant that the decline rates of gas production for different scenarios were almost the same. After the time, t>100 days, the differences become smaller with increase in time. In addition, it was also observed that with decrease in formation thickness, the area enclosed by the curve became smaller, which meant the cumulative production reduced.Figure 5
Effect of formation thickness on production decline curves.From Figure6, it was noted that permeability affects the whole period of gas production. The higher the fracture permeability, the higher the initial production rate. However, for a reservoir with higher fracture permeability, production rate would decrease drastically, and during the middle-late stages the production rate would be exceeded by a reservoir with lower fracture permeability. Moreover, a well with low fracture permeability has a long stable production period.Figure 6
Effect of fracture permeability on production decline curves.Figure7 shows the effects of reservoir radius on gas production performance of HDWs. The difference of gas production at late stage was increasingly wider. It was noted that with increase in gas reservoir radius, the production curves declined more gently, which meant a gas reservoir with bigger radius would have a longer stable production period.Figure 7
Effect of gas reservoir radius on production decline curves.Figure8 shows the effects of interporosity flow coefficient between fractures and vugs on gas production performance of the HDW. It can be seen that the coefficient primarily affected early stage of gas production. For the time, t<10 days, with the interporosity flow coefficient between fractures and vugs smaller, the initial production also became smaller and the production rate decreased relatively slower. For the time, t>10 days, there was no difference between the production rates of the three scenarios. From Figure 9, it was observed that vugs storativity ratio affected middle-late stage of gas production. During 10<t<10000 days, with vugs storativity ratio smaller, production rate decreased more rapidly. However, after t>10000 days, the difference in the production curves was contrary to the previous one. Meanwhile, it was obvious that the extent of production variation is less drastic than the middle stage.Figure 8
Effect of interporosity flow coefficient between fractures and vugs on production decline curves.Figure 9
Effect of vugs storativity ratio on production decline curves.Figure10 shows the effects of the angle of inclination on gas production performance of the HDW. It can be seen that the angle of inclination also affects the entire stage of gas production. Due to the effects of different inclination angle on contact area of wellbore, the cumulative gas production increases with increase in inclination angle.Figure 10
Effect of inclination angle on production decline curves.According to Figures5–10, the production decline curve was affected by many parameters and the intensity of production variation was different at different stages of gas production history. By extracting production rate data at different times (t=10days, 100days, 1000days) for adjacent parameters in Table 3, production difference percentage (Δqper) could be obtained which could be helpful in comparing the sensitivity of these parameters.Table 3
Effect of parameters att=10 days, 100 days, 1000 days on the difference in production rate percentage.
parameters  10 days 100 days 1000 days Δ q p e r, % Δ q p e r, % Δ q p e r, % k f, md 1/0.5 22.1 18.4 18.2 r e, m 800/400 3.8 12.8 20.9 λ c 0.0002/0.00002 4.6 1.1 1.2 ω c 0.1/0.01 3.7 7.5 5.9 θ, degree 80/40 23.1 21.8 18.7Based on Table3, trends of difference in percentage of parameters are shown in Figure 11. It can be seen that the bigger the difference in percentage value, the more sensitive the parameters. Inclination angle (θ) and fracture permeability (kf) seemed to be the two most sensitive parameters; reservoir radius (re) was the third, vugs storativity ratio (ωc) and interporosity flow coefficient between fractures and vugs (λc) were the least two sensitive parameters.Figure 11
Trends of percentage difference of parameters by extracting production rate att=10 days, 100 days, 1000 days.


## 5. Conclusions

This research investigated the production decline behavior of highly deviated well (HDW) in naturally fractured-vuggy carbonate gas reservoir. The validity of the proposed model was verified by matching it with the production history of a field case from the Arum River carbonate gas reservoir in Turkmenistan. Furthermore, the effects of relevant parameters on production performance were studied. The main conclusions drawn are as follows:(1)
The proposed production semianalytical model could be used to predict the gas production trend with a constant bottom-hole pressure.(2)
Pressure type curves for HDW in naturally fractured-vuggy carbonate gas reservoir were divided into five flow stages. Among them, the second stage was dominated by the inclination angle (θ) of the HDW. The third stage was dominated by interporosity flow between fractures and vugs. The fourth stage was dominated by interporosity flow between fractures and matrix. However, the third and fourth stages could interfere with each other, depending on the values of λc and λm.(3)
Using the proposed model, sensitivity analyses were conducted and it was found that formation thickness (h), fracture permeability (kf), inclination angle of the HDW affected the entire stage of gas production history. Interporosity flow coefficient between fractures and vugs (λc) affected the early stage of gas production history. Vugs storativity ratio affected middle-late stage of gas production history. Gas reservoir radius affected late stage of gas production history. Moreover, inclination angle (θ) seemed to be the most sensitive parameter, fracture permeability (kf) was the second most sensitive parameter, followed by reservoir radius (re), while vugs storativity ratio (ωc) and interporosity flow coefficient between fractures and vugs (λc) were the last two sensitive parameters.​



---

*Source: 6919176-2019-07-29.xml*