Abstract
Fracture-vuggy oil reservoirs contain micrometer-sized dissolved secondary pores, fractures, and meter-sized karst caves. Generally, the researches on the flow mechanism of fractured-vuggy reservoirs are based on the assumption of karst cave-fracture-matrix triple medium that gives a partial differential equation which the pressure satisfies. This triple medium flow equation lays the basics for the well testing theory of fractured-vuggy reservoirs. However, the triple medium equation cannot reflect the actual existence of the karst cave volume in geology. Recently, we proposed a theory of fracture-cave seepage flow coupled with vug pressure wave which takes the vug size into consideration and thus can better describe the fracture-vuggy reservoirs. On the basis of this theory, this work combines the change of formation pressure with time, Duhamel’s theorem, and well testing theory to obtain the mathematical expression of time-varied bottom hole pressure and flow rate in Laplace space. Then, the dynamic inflow performance relationship (IPR) curve of fractured-vuggy reservoir is obtained by inverse Laplace transformations. This dynamic IPR curve can predict the productivity for fractured-vuggy reservoirs using parameters from well test interpretations of pressure recovery. The correctness of this productivity prediction method is verified by filed fractured-vuggy oil well. This productivity prediction method not only utilizes and expands the application of well test interpretation data but also shortens the test time and reduces test costs, which is important for oil development.
1. Introduction
Fractured-vuggy carbonate reservoirs contain micrometer-sized dissolved secondary pores, fractures, and meter-sized karst caves (the large vugs). Fractures and caves are the main storage space for the oil resources [1]. Pores-fractures-caves make the development of fractured-vuggy reservoirs accompanied by strong heterogeneities and discontinuities, which is different from the sandstone reservoirs [2, 3]. Productivity prediction refers to the technical evaluations of the daily or yearly productivity of oil wells. More specifically, well productivity refers to the relationship between bottom hole pressure and surface production during single well production. It is a core parameter in oil and gas development. The prediction results will affect the following working system, drilling workload, investment scale, etc., which indicates that the productivity prediction is very important to the successful development of oil reservoirs.
Many researches have been carried out on productivity well testing methods relating to conventional reservoirs [4, 5]. Backpressure well testing, isochronal well testing, and modified isochronal well testing methods lay the basic theories for productivity testing, and subsequent stable well testing studies are based on the application of these theories. The backpressure well test method developed by Rawlins and Schellhardt requires that both production and bottom hole flow pressure achieve to be stable with time. Therefore, the time cost of this method is relatively long (usually in several or ten days above) and the amount of released gas into air is quite large, which decreases profit to the producers [6]. The isochronal well testing method developed by Cullender and Smith does not require the production and pressure to be stable, which not only reduces the testing time but also reduces the amount of air venting [7]. However, this method requires multiple well shut-in manipulations, and each time the well is shut-in, it must be stabilized and restored to the original formation pressure. The operation procedure is more complicated than the backpressure well test and requires more time. The modified isochronal well test method does not require the original formation pressure restored when shutting in the well. This method is especially suitable for low-permeability formations which needs more time to restore original formation pressure [8–10]. The one-point method has better economic benefits for gas well productivity testing [11]. The open flow capacity can be quickly obtained once the formation pressure, the bottom hole pressure, and the corresponding gas output are measured under one working cycle when using this test method. However, these productivity prediction methods assume constant pressure, flow rate, or time, which actual oil wells have synchronous changes of both pressure, flow rate, and time.
Except for the productivity well testing methods, many theoretical deductions of productivity equation have emerged due to the high cost of well testing which is based on the theory of stable seepage. Vogel used a computer to simulate the saturated oil flow processes and concluded the relationship between bottom hole pressure and flow rate [12]. Vogel’s method is suitable for oil well productivity prediction in two-phase oil and gas flow. For high-yield oil and gas wells, the turbulent flow has nonnegligible impact. Fetkovich proposed three different forms of equations that can better describe the inflow performance of such wells [13]. Wiggins obtained the dimensionless oil equation and water equation by linear regression on the results of the simulator and used them to describe the oil-inflow dynamic characteristics of the oil well during the three-phase flow of oil, gas, and water [14].
For fracture-vuggy reservoirs, the triple medium model is adopted in many researches but there are seldom researches on the fracture-vuggy reservoir’s productivity prediction. Camacho-Velazquez et al. proposed the triple-porosity/dual-permeability model to include the interactions between matrix, fractures, and vugs [15]. This three-continuum further developed to multiple-continuum models [16]. Huang et al. proposed a numerical model for immiscible two-phase flow in fractured karst reservoirs based on equivalent continuum representation [17]. Wu and Ge proposed the three-porosity model of naturally fractured reservoirs in both infinite and finite systems by using the Laplace transform and the finite Hankel transform [18]. Xu et al. carried out the dual medium mathematical productivity prediction model of fracture cave carbonate reservoir in buried hill and analyzed the influencing factors of productivity [19]. Wang et al. used the triple medium model and the Laplace transform to obtain the productivity equation of a single well in a fracture-vuggy reservoir [20]. This triple medium model treats the cave as a porous medium but not as a void vuggy, then uses the interporosity flow equation to describe the matrix-fracture-cave flows between them [16]. However, this model ignores the flow in the vuggy and simplifies it to cross flow which is not consistent with the actual production. Du et al. proposed an analytical model for fractured-vuggy carbonate reservoirs considering the vuggy volume using the wave propagation theorem [21]. The wellbore pressure history calculated by Du et al.’s model matched well with the field measured data from Tahe reservoirs in China.
Since many wells belong to depletion development since no extra gas or oil beyond the reservoir is supplied to keep the formation pressure from going down [22], Arps proposed three types of productivity decline curves with time: the exponential, harmonic, and hyperbolic declines through summarizing a large amount of field data [23]. Everdingen and Hurst, Fetkovich, and Blasingame et al. further enriched and developed the productivity decline curve analysis from the seepage equation and finally developed them as the production data analysis method [24–26]. These methods deduces several formation parameters through data analysis of the everyday production pressure and flow rate change. Moreover, the productivity testing and productivity equations are based on steady seepage theory and do not consider the impact of time on productivity. The production decline curve considers the time effect, but it is based on the assumption of constant bottom hole pressure. Jiang and Lu et al. proposed the concept of dynamic IPR curve between well production and bottom hole pressure at different time by constantly changing the time to calculate the productivity using the interpretation results of pressure recovery well test or production data [27, 28]. The dynamic IPR curve considers the relationship between pressure and flow rate with time, but it assumes that the formation pressure remains unchanged.
Therefore, currently, there are seldom appropriate methods to predict the productivity of fractured-vuggy carbonate reservoir wells because actual reservoirs usually have synchronous changes of both pressure, flow rate, and time [29, 30]. The seepage flow of fractured-vuggy reservoirs is different from that of conventional reservoirs because of the large and rapid drainage of the vuggy cave in the fractured-vuggy formations. There is no special method for predicting productivity for fractured-vuggy wells. Therefore, this work proposed a new predicting productivity method for fractured-vuggy reservoirs considers the mathematical well flow equations, the time change of bottom pressures, the flow rate, and the averaged formation pressure (i.e., the stratum static pressure) caused by long-term production. Field fractured-vuggy oil wells are taken as an example for the application of this productivity predicting method and show reasonable accuracy.
2. Theory and Methodology
In this section, firstly, we introduce the mathematical model for fractured-vuggy reservoirs which is a different traditional triple medium model that ignores the dynamic flow in the vugs; secondly, we introduce the theory of getting time-varied average formation pressure due to the production that causes the formation pressure decreasing; thirdly, we introduce the new method to get the productivity predictions for fractured-vuggy reservoirs with time-varied flow rates or bottom hole flow pressures.
2.1. Mathematical Model for Fractured-Vuggy Reservoirs
Fractured-vuggy reservoirs are spatially discrete media, characterized by complex storage space, coexistence of fractures and large-scale karst caves, complex spatial configuration of karst fractures, and diverse filling types and filling levels [31, 32]. In our previous work [21, 33–35], we proposed a flow model that coupled pressure wave conduction in the cave and seepage flow and defined the wave coefficient, shape factor, and damping coefficient of the cave. This model as shown in Figure 1 has the following assumptions [35]: (1) the shape of the vug connected with the wellbore is a cylinder, and the wellbore sits in the top center of the vug; (2) the flow from the vug to the well is only vertical. The formation flow obeys Darcy’s law and is isothermal; (3) the formation is isotropic and cylindrical. The reservoir contains only oil; (4) the outer formation is composed of a matrix; the formation permeability and porosity and the fluid compressibility and viscosity are constants.

To build the mathematical equations for the fractured-vuggy reservoir, these dimensionless variables are firstly defined:
is the dimensionless pressure for formation , , well bottom hole, and vugs; is the dimensionless time; is the dimensionless distance; is the dimensionless height weighted; is the dimensionless wellbore storage constant; is the dimensionless vug storage constant; is the dimensionless cave damping coefficient; is the dimensionless cave fluctuation coefficient; is the dimensionless friction coefficient; is the dimensionless fluctuation coefficient.
The definitions and units of these variables are as follows: is the original formation pressure (MPa); are the pressures in formation , , well bottom hole, and vugs, respectively; is the viscosity of the fluid (mPa·s); and are the wellbore and cave storage constants (m3/MPa), respectively; and are the skin factors of the wellbore and cave; is the daily production (m3/day); is the production time (day); is the fluid volume factor; and are the wellbore and cave radius (m), respectively; is the fluid density (kg/m3); is the formation permeability (μm2); is the effective thickness of the formation (m); is the initial velocity (m/s).
Using the dimensionless variables defined above, the following dimensionless equations [35] using Laplace transform are obtained as
The general solution of in Laplace space for a closed circular reservoir can be obtained as where
where and are the first kind modified Bessel functions with order 0 and order 1, respectively; is the second kind modified Bessel functions with order 1. is the function, and is the Laplace variable.
The solution of in real space can be obtained using the above Laplace solutions by the Stehfest numerical inversion method.
In this equation, is an even number that generally ranges from 8 to 16.
2.2. Time-Varied Average Formation Pressure
According to the principle of material balance, for a sealing oil formation, if there is no external energy supplement, the formation pressure will decrease with continuous production. The averaged formation pressure is obtained by integrating the pressure distribution over the area inside the formations as shown in Figure 2. The pressure distribution is usually represented by the pressure funnel curve as shown in Figure 3.


For a given production time, the formation pressure distribution can be obtained by solving the formation seepage mechanic equations. The average pressure can be defined as where is the time-varied average formation pressure, the production time, is the drainage area, and is the pressure distribution at any location point in the formation at time .
If the reservoir formation is circular, using circular coordinates, equation (5) can be simplified to where is the pressure distribution at any point in the formation at time , is the formation radius, and is the well radius.
The average formation pressure model assumes that (1) the formation and fluid have microcompressibility, (2) the gravity and skin factor and well storage effect are neglected (since only the average formation pressure is considered), (3) there is homogenous circular formation with one well in the center of it, and (4) the outer is closed boundary with radius .
The dimensionless bottom hole pressure during the average formation pressure calculation is defined as
The dimensionless production time is defined as
The dimensionless radius and are defined as
Apply the Laplace transform of mathematical equation for circular closed formation ( is the Laplace variable), we get
The above equation has a solution of as where and are the first kind modified Bessel functions with order 0 and order 1, respectively; and are the second kind modified Bessel functions with order 0 and order 1, respectively.
The process of inversion Laplace transform of is quite complex [36], but it has analytic solutions for in real space, and using the integral theory on the mathematical pressure expressions [37], the average formation pressure is calculated as where is the decreased average formation pressure along with well production, is the dimensionless production time, is the original formation pressure, is a small quantity that characterizes the radius of the well, is the production drainage area, is the volume factor of fluid, is the formation permeability, is the comprehensive compression factor, is the formation porosity, is the effective stratum thickness, and is the daily production.
Substituting the dimensionless time into equation (12) and considering the total volume produced in general form of , then equation (12) can be simplified to the expression of formation average pressure versus produced volume: where is the total produced fluid volume under standard ground conditions and is a time-varied parameter.
Equation (13) indicates that the time-varied average formation pressure does not include flow parameters such as permeability or fluid viscosity, which indicate that the average formation pressure has no relationship with the seepage process of formation fluids. This important feature can be further applied to the fracture-vuggy reservoirs to calculate its average formation pressure: firstly, the mathematical model of the fracture-vuggy reservoirs makes it difficult to obtain the analytical solution of average formation pressure by the integral theory because the pressure distribution can only be obtained by the Stehfest numerical method; secondly, this equation is derived from one well without the vug, which leads us to consider whether the pressure in the vug could influence the average formation pressure; thirdly, the pressure transmits fast from the vug to the well in the presence of waves and equilibrates the heterogeneous pressure inside the vug, making the vug pressure be treated as one point and neglected. Therefore, we can use equation (13) to recalculate the average formation pressure after production for fracture-vuggy reservoirs.
2.3. Productivity Prediction Model for Fractured-Vuggy Reservoir
Duhamel’s theorem is an important method in the analytical solution of partial differential equations [38]. It reduces the problem of boundary conditions and nonhomogeneous terms that change with time to the problem of partial differential equations and then solves them, thus simplifying the solution. For oil and gas development, daily production and pressure always change with time. According to Duhamel’s theorem, we firstly solve the pressure solution for one unit flow rate (represented by subscript ) and then obtain the real pressure solution through the integral of convolution between and the real flow rate . The general solution expression of Duhamel’s theorem is where is the bottom hole pressure, is the original formation pressure, and is the pressure difference per unit flow rate. Considering that the average formation pressure changes during the long-term production of the formation, in equation (14) should be replaced by the average pressure in equation (13), and the cumulative volume produced can be expressed by flow integral, so that equation (14) becomes where is defined as pressure decreasing coefficient (MPa/m3).
Equation (15) is a general form because it is directly expressed by the partial differential equation and the average pressure expression. This expression considers both the influences of simultaneous changes of flow rate and formation pressure on bottom hole pressure.
Because the mathematical solution for fractured-vuggy reservoirs can only be obtained by the numerical inversion of the Laplace transform, then equation (15) in real space cannot be used directly to get the productivity, i.e., the relationship between pressure and flow rate. We have to deal with it from the Laplace space. According to the Laplace properties of integral and convolution, the bottom hole pressure difference in the Laplace space from equation (15)can be expressed as where is the bottom hole pressure difference, the overbar indicates of the Laplace transform expression where is the Laplace operator, is the function in real space, and is the corresponding function in the Laplace space. acts on , , and , respectively.
When the flow rate is one unit, according to the definition of dimensionless pressure, the bottom hole pressure can be expressed as
Perform the Laplace transform of equation (17), the unit bottom hole pressure in the Laplace space is
In Section 2.1, we have obtained the expression of in equation (2). Therefore, substituting equation (18) and equation (2) into equation (16), the relationship between bottom hole pressure of fractured-vuggy reservoirs under time-varied flow rate is obtained as follows:
This is the Laplace equation of productivity prediction for fractured-vuggy reservoirs. This new equation (19) is a step closer to actual production by considering both the time-varied bottom hole and formation pressures and flow rate. Since there are no analytical solutions for the pressure calculations of fractured-vug reservoir, the detailed calculation of productivity for fractured-vug reservoirs is obtained by the Stehfest numerical method.
3. Applications
Using the method in this work, equation (19) shows that the calculation of productivity requires parameters such as permeability and cave volume as known. These parameters can be obtained in three ways relating to experiments or tests. Each experiment or test has its own meaning and effects of productivity prediction. (1) Obtain the data during the exploration period and from the drilling-log information, which can estimate the future development potential and determine the development value of the exploration area; (2) obtain these parameters for productivity calculation through the interpretation of the well test analysis data in the exploration and test stage. The oil development plan can be determined more accurately from the production capacity calculation; (3) obtain these parameters through the monitored pressure and flow rate during production, these parameters are obtained through production data analysis technology to predict the production capacity, and the development system can be optimized based on the production capacity data. This work will use the second way, i.e., the reservoir well test pressure recovery data to predict productivity and compare the predicted results with actual test data to verify the correctness of our method.
According to the pressure recovery field test data of one fractured-vuggy well in Xinjiang, China, the basic parameters of formation, well, and fluid are shown in Table 1 [35]. The measured historical flow data for the same oil well which is obtained from Sinopec is shown in Figure 4.

Using the mathematical model of fractured-vuggy reservoir, the interpretation results of this pressure recovery field data are shown in Table 2. Through this, then the IPR can be calculated. Equation (19) gives the relationship between flow rate and bottom hole pressure in the Laplace space. Using the Laplace inverse transformation, the relationship between flow rate, bottom hole pressure, and time can be obtained. If the planned production duration (days) is given, then the relationship between flow rate and bottom hole pressure, i.e., the dynamic IPR curve can be calculated and drawn. Here, since the flow rate-pressure change is related to time, so it is called the dynamic IPR curve, compared with traditional IPR curve that will not change with time [27, 28].
According to the data provided in Tables 1 and 2, the dynamic IPR curves with continuous production time for 1 year, 2 years, and 3 years are calculated, respectively. Figure 5 and its enlarged part view Figure 6 show the calculation results.


The measured historical flow chart of one fractured-vuggy well in Xinjiang as shown in Figure 4 indicates that this well’s production within 3 years ranges mainly from 100 to 130 m3/d, while the bottom hole flow pressure lowered to be 83 MPa. From the calculated dynamic IPR curves in Figure 6: if we want continuous constant production for 1 year and want the bottom hole flow pressure no less than 83 MPa, then the flow rate should be no higher than 135 m3/d. Similarly, if we want to produce for 2 years, 3 years, and the bottom hole flow pressure no less than 83 MPa, then the flow rate should be no higher than 120 m3/d, 110 m3/d, respectively. Therefore, the predicted production matches reasonably well with the measured production and shows that this proposed productivity prediction method is effective.
Considering that these oil wells in Xinjiang have no extra energy supplement, the formation static pressure will decrease with production. According to the measured static pressure data of this oil well, the decrease rate of the formation static pressure is about 2.4 MPa/year. Here, we calculated the dynamic IPR curves under variable static pressure with continuous production for one year as shown in Figures 7 and 8 (partially enlarged view of Figure 7).


The measured formation static pressure for this oil well is 86.8 MPa for the 1st year as shown in Table 2, and the measured decrease rate of the formation static pressure is about 2.4 MPa/year, i.e., the formation static pressure for the 2nd year is 84.4 MPa, and 82 MPa for the 3rd year. Then, the flow rate under constant pressure difference of 2.8 MPa (bottom hole flow pressure at 84 MPa from Figure 8) for one year is calculated to be 100 m3/d. The dynamic IPR curves are parallel to each other at a low flow rate indicating that the sufficient energy is maintained due to high static pressure. As shown by the dotted lines in Figure 8, when static pressure in the 2nd year is 84.4 MPa, a pressure difference of 2.8 MPa could still produce a 100 m3/d flow rate. The same analysis is applicable for the 3rd year with a static pressure of 82 MPa, and the predicted flow rate at pressure difference of 2.8 MPa is 100 m3/d.
The production flow rate over time when this oil well is produced at 1st year at constant pressure difference of 2.8 MPa is shown in Figure 9. The predicted flow rate of this oil well for one year is around 110 m3/d and is consistent with the measured results as shown in Figure 4, i.e., when the production time is less than 400 days, the actual production is around 110 m3/d. This indicates that our productivity prediction method for fractured-vuggy reservoirs agrees well with the actual field production data.

4. Conclusions
In this work, we combined the time-varied average formation pressure, Duhamel’s theorem, and fractured-vuggy reservoir model to develop a new equation and method of productivity predictions. The following conclusions are drawn through this research: (1)Based on the coupling theory of vug pressure wave and formation pressure seepage, a new method was proposed to get productivity predictions for fractured-vuggy reservoirs with time-varied flow rates or bottom hole flow pressures. The application situations of this method are closer to the real field production than traditional methods(2)The fractured-vuggy reservoir model, average formation pressure theory, and Duhamel’s theorem are implemented together to get the flow-pressure relationship (the IPR curves) under different production times(3)The proposed method is applied and verified by the productivity prediction of one fractured-vuggy oil well in Xinjiang, China. The predicted flow rate is in reasonable agreement with the actual measured data. Since only the formation parameters which be obtained through one-time shut-in pressure recovery operation are needed, this method avoids the frequent operations of well open and shut-in for traditional productivity testing and therefore saves both cost and time
Considering the unstable production plans in oilfield, this work provides a new way of productivity prediction for fractured-vuggy reservoir and could help on the optimization of fractured-vuggy reservoir development plans.
Data Availability
The historical flow data of one fractured-vuggy oil well used to support the findings of this research is available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by “Research on quality improvement and speed up drilling and completion technology of No. 5 Fault Zone in ShunBei Area” with contract No. P20002 held by Sinopec scientific research project in 2020.