#### 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.

#### 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

Figure 1 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 to .(6)The HDW produces at a constant rate or at constant pressure .(7)The permeability in horizontal direction and in vertical direction are different.

#### 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 Appendix A), the governing equation for the triple-porosity system in radial cylindrical coordinate is given as follows:

For the matrix system the equation is

For the vugs system the equation is

where , , is the pseudopressure of fractures system, matrix system, vugs system, respectively, ; , , is porosity of fracture system, matrix system, vugs system, respectively, fraction; , , is permeability of fracture system, matrix system, vugs system, respectively,* md*; is vertical fracture permeability,* md*; , , is compressibility of fracture system, matrix system, vugs system, respectively, ; is gas viscosity, ; , are shape factors of matrix and vugs, ; and is pseudotime, day.

The initial pseudopressure is assumed uniform for the triple-porosity system:

The outer boundary is closed:

At the inner boundary, there is a continuous point source () with a constant gas rate :

where is fracture horizontal permeability, .

At the top and bottom, the boundaries are impermeable:

where is formation thickness, .

Based on the definitions of dimensionless variables in Table 1, the governing equations and the initial and boundary conditions (Eq. (1)-(7)) can be transformed into dimensionless terms (Eq. (8)-(14)).

The dimensionless governing equation of fracture, matrix, and vugs systems are

Transformed initial condition:

Transformed outer boundary condition:

Transformed inner boundary condition:

Transformed top and bottom boundaries:

##### 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).

where

where is modified Bessel function of zero order of the second kind; is modified Bessel function of order one of the second kind; is modified Bessel function of order zero of the first kind; 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 the , , -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:

where

Based on the principle of superposition, the pressure solution with wellbore storage and skin effect in Laplace domain can be obtained as

where* C* is wellbore storage coefficient, ; 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:

Equation (25) can then be inverted to obtain 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 Table 2. These parameters are obtained by various methods in the field. In particular, it should be noted that formation radius (), fracture permeability (), interporosity flow coefficients between different media systems (), and initial reservoir pressure () can be obtained through interpretation of well testing data. And the porosity of three different media () can be obtained by combining well logging data with core CT scanning.

Figure 3 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 () between the remaining actual data and model data is only 4.19%. The average relative error () is evaluated by

where is simulated production rate from the proposed model, ; is field production rate, .

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. Figure 4 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 and . The fifth stage was dominated by closed circle boundary, which manifested on the derivative curve with a slope of 1.

##### 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.

Figure 5 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, 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, 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.

From Figure 6, 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 7 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 8 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, 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, 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 days, with vugs storativity ratio smaller, production rate decreased more rapidly. However, after 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 10 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.

According to Figures 5–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 () could be obtained which could be helpful in comparing the sensitivity of these parameters.

Based on Table 3, 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 () seemed to be the two most sensitive parameters; reservoir radius () was the third, vugs storativity ratio () and interporosity flow coefficient between fractures and vugs () were the least two sensitive parameters.

#### 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 and .(3)Using the proposed model, sensitivity analyses were conducted and it was found that formation thickness (), fracture permeability (), inclination angle of the HDW affected the entire stage of gas production history. Interporosity flow coefficient between fractures and vugs () 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 () was the second most sensitive parameter, followed by reservoir radius (), while vugs storativity ratio () and interporosity flow coefficient between fractures and vugs () were the last two sensitive parameters.

#### Appendix

#### A. Calculation of Pseudopressure and Pseudotime

In (1)- (7), , , is the pseudopressure of the fracture system, matrix system, and vugs system, respectively, ; is pseudotime, day. The pseudopressure and pseudotime can be given as follows:

#### B. Calculation of Continuous Point Source Solution

Equations (8)-(14) can be transformed into Laplace domain by the Laplace transformation.

The dimensionless governing equations of fractures, matrix, and vugs systems in the Laplace domain are as follows:

The outer boundary condition is transformed as

The inner boundary condition is transformed as

The top and bottom boundaries are transformed as

where

To eliminate the variable in the governing equation, (B.1)-(B.5) can be transformed by Fourier cosine transform. The Fourier cosine transform and inverse Fourier cosine transform are given as follows:

where

By employing the Fourier cosine transform (Eq. (B.8)), Eqs. (B.1)-(B.5) can be transformed as follows:

where

The outer boundary condition is transformed as

The inner boundary condition is transformed as

Then (B.11) can be transformed to a modified Bessel function of zero order, like the following:

The general solution of (B.16) is given as follows:

Based on the outer and inner boundary conditions (Eq. (B.14)-(B.15)), the solution in Laplace domain can be obtained by inverse Fourier cosine transform as follows:

#### Nomenclature

Wellbore storage coefficient, | |

: | Total compressibility, |

Average relative error, % | |

Formation thickness, | |

Permeability of vugs system, | |

: | High deviated well length, |

Pseudopressure, | |

: | Well bottom-hole pseudopressure, |

Pressure, | |

: | Well bottom-hole pressure, |

: | Pressure at standard condition, |

: | Production rate from point source, |

: | Gas production rate, |

: | Simulated production rate from the proposed model, |

: | Field production rate, |

Radial distance, | |

: | Wellbore radius, |

: | Formation radius, |

Skin factor | |

Laplace transform variable | |

Pseudotime, day | |

: | Pseudotime, day |

Reservoir temperature, | |

: | Temperature at standard condition, |

Directional coordinates | |

: | Distance of mid-perforation in x, y, and z coordinates, |

Z-factor of gas, dimensionless | |

: | Shape factors of vugs and matrix, |

Interporosity flow coefficient, dimensionless | |

Storativity ratio, dimensionless | |

Inclination angle, degree | |

Porosity, fraction | |

: | Gas viscosity, mPa·s |

Constant, |

*Subscript*

Vugs system | |

Fractures system | |

Vugs system | |

Horizontal direction | |

Initial condition | |

Initial condition | |

Vertical direction | |

Dimensionless. |

*Superscript*

Laplace domain | |

Fourier domain | |

Field production data. |

#### Data Availability

The field data from the Arum River Basin in Turkmenistan used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors would like to acknowledge the support provided by Xin Zhao and Jun Shi in the initial stages of this study. This article is funded by National Science and Technology Major Project of China (*Grant Nos. *2017ZX05030-003* and *2017ZX05009-005).