A Practical Methodology for Production Data Analysis in Carbonate Reservoirs Using New Decline Type Curves

Carbonate reservoirs typically have complex pore structures, so the production wells typically have high production in the early production stage, but they decline rapidly. It is highly challenging to achieve accurate interpretation results. In this paper, a new and practical methodology for production data analysis of fractured and fractured-vuggy carbonate reservoirs is proposed. Firstly, analytical solutions to characterize the different multipore media and simulate transient production behavior of fractured and fractured-vuggy carbonate reservoirs during the transient flow regime are presented. Then, a new function f ðqDÞ and f ′ðqDÞ that related to the dimensionless production rate is introduced, and a series of new decline type curves are drawn to make a clear observation of different flow regimes. In addition, the effects of the storativity ratio, interporosity flow coefficient, skin factor, and dimensionless radial distance of external boundary on production performance are also analyzed. Finally, two example wells from the fractured and fractured-vuggy carbonate reservoirs are used to perform rate decline analysis with both the Blasingame type curves and the new type curves. The validation of the new method is demonstrated in comparison to the results of well test interpretation. The results show that the curves of 1/f ′ðqDÞ vs. tD are ∧-shaped for dual-porosity reservoirs and M-shaped for triple porosity reservoirs and also indicate that the interpreted parameters such as permeability, skin factor, storativity ratio, and interporosity flow coefficient using new decline type curves are aligned well test interpretation. In correlation with other traditional well test analysis, this approach effectively reduces the multisolution probability of interpretation.


Introduction
Carbonate reservoirs, which approximately account for 60% of the world oil reservoirs, represent a significant amount of oil and gas reserves and have a great potential to be exploited [1]. With further exploration and development, these fractured (see Figure 1(a)) and fractured-vuggy (see Figure 1(b)) carbonate reservoirs have been attracting considerable attention in recent years [2][3][4][5]. The complex structures of the reservoir are characterized by the presence of different porous media and transmissibility [6]. Production wells of carbonate reservoirs typically have high production in the early production stage but they decline rapidly, leading to the difficulties of further exploitation [7][8][9][10][11]. At present, depletion development is mainly applied to extract oil and gas from carbonate reservoirs, and the pressure drawdown is quite large to have a considerable production rate. Nevertheless, it is difficult to predict the production and analyze production data for fractured and fractured-vuggy carbonate oil reservoirs.
In order to characterize the complex pore structures and flow features of fractured and fractured-vuggy carbonate oil reservoirs, the concept of multiple-porosity is proposed. For the dual-porosity, Barenblatt et al. [12] hypothesized that there are two continuous media, matrix and fracture. The porosity of fracture is much less than that of matrix, but the permeability of fracture is much greater than that of matrix. Warren and Root [13] proposed the idealized dual-porosity model and solved the governing equation of pseudosteady flow analytically. They assumed that the fracture system evenly divides the matrix rock. The matrix system which is in a hexahedral shape is homogeneous and isotropic. Kazemi [14] proposed a new model which is assumed that the reservoir is composed of horizontal fractures and matrix, and the fracture system evenly divides the matrix system. De Swaan [15] developed a model and assumed that the configuration relationship between the matrix system and the fracture system is similar to that in the Warren-Root model, except that the matrix rock blocks are round sphere. For the tripleporosity, Abdassah and Ershaghi [16] firstly proposed the triple-porosity model, which considers that the reservoir is composed of three media: matrix, fracture, and vug. The model is assumed that the fracture communicates directly with the wellbore, and the matrix and vug simultaneously have an unsteady flow to the fracture. Jalali and Ershaghi [17] considered the heterogeneity of the reservoir and analyzed the unsteady pressure behavior of such reservoirs. Al-Ghamdi and Ershaghi [18] developed the triple-porosity/dual-permeability model, which proposed two types of flow units to differentiate between microfractures and macrofractures. Based on the theories and models of these pioneers, indepth research of the dual-porosity and triple-porosity reservoirs can be found in the recent literature [19,20].
In recent years, empirical, numerical, and analytical methods have been widely used for production data analysis of fractured and fractured-vuggy carbonate reservoirs [20][21][22][23][24]. Because of its relative simplicity, empirical methods are commonly applied in the actual oil field to conduct history matching and make forecasts for production wells [21,25,26]. Arps [25] decline curve analysis is a classic empirical method used mainly to predict recoverable reserves and future production rates. However, this method is based on empirical observations of production decline. Furthermore, the decline rate greatly dominates the decline curve analysis. Conventional methods have been the use of numerical simulation for performing production data analysis of fractured and fractured-vuggy carbonate oil wells [22,27,30,31]. However, it is time-consuming in gridding and simulating, which makes numerical models less attractive. Analytical models not only capture the reservoir architecture of carbonate reservoirs but also have much higher computational efficiency than numerical models, which have recently attracted more and more attention [25][26][27][28][29][30][31]. Fetkovich [32] expanded Arps's approach to cover the transient flow condition using analytical flow equations. However, this method is only suitable for a well producing at constant bottom-hole pressure (BHP) condition while the Pressure-Volume-Temperature (PVT) properties are not considered. Blasingame et al. [33] proposed a rate decline analysis method that accounts for variable BHP and production rate conditions. This allowed the analysis of well production rate with variable BHP and variable production rate, by reducing them to an equivalent constant-rate data. In addition, the PVT properties are also incorporated in the models. Therefore, Blasingame type curves have been widely used for field data analysis. Agarwal et al. [34] expanded Fetkovich's and Blasingame's work to develop a new type decline curve based on well test definitions, as opposed to the dimensionless variables in the Blasingame method. By such, an approach to the rate decline analysis is to reduce the multisolution to certain content. In practice, the boundary-dominated regime should be observed on the type curve to avoid multisolution. Based on these theories and methods, extensive researches and field applications have been carried in the literature [35][36][37][38][39]. However, when the wells exhibit a long transient flow regime and the boundary-dominated regime could not be observed, especially the flow behavior characteristics of fractured and fractured-vuggy carbonate reservoirs, these methods are generally led to multisolution. Therefore, proposing an effective method to capture the production performance of wells during the transient flow regime is of significant importance [40].
Firstly, a novel methodology is presented to analyze production data of wells from fractured and fractured-vuggy carbonate reservoirs during the transient flow regime. Then, a traditional dual-porosity and triple-porosity model is applied to model the fractured and fractured-vuggy carbonate reservoirs, respectively. With Laplace transform, the analytical solution of the mathematical model is obtained. In order to develop the new type curves (see Table 1), functions related to the dimensionless production rate are introduced. Furthermore, a series of new type curves are drawn to make a clear observation of different flow regimes. The effects of some critical parameters on production performance are also analyzed with the proposed model. In addition, the new method is benchmarked against the well test interpretation. Finally, two field cases are used to demonstrate the practicability of the method. The main contribution of this paper is the provision of a novel approach for production data analysis of fractured and fractured-vuggy carbonate oil wells during the transient flow regime. Compared to other methods, the accuracy of the proposed method for production data analysis is higher. This new method provides another choice for field engineers and related research and further reduces the multisolution probability by comparing among the different interpretation methods. Although the traditional multiple media model has been studied in this paper, this new method is also applicable to other complex fractured  2 Geofluids reservoirs. The corresponding research work will be carried out in the next step.

Methodology
2.1. Physical Model and Assumptions. It can be seen from Figure 2(a) that the fractured carbonate reservoir is modeled with a dual-porosity model, which is structured by fracture and matrix systems. The single-phase oil flows from the matrix system to the fracture system, and then to the wellbore (see Figure 2(b)). The fractured-vuggy carbonate reservoir is modeled with a triple-porosity model, which is structured by fracture, vug, and matrix systems (see Figure 2(c)). The oil firstly flows from the vug system to the fracture system, then from the matrix system to the fracture system, eventually to the wellbore (see Figure 2(d)). The matrix-fracture interporosity flow and the vug-fracture interporosity flow in the reservoirs are described by the pseudosteady state model [9,16,19]. The radial cylindrical dualporosity and triple-porosity reservoirs are considered in which a single production well is located at the center, completely penetrating the formation. Some simplifying physical model assumptions for the derivation of the governing equation are as follows: (1) Total compressibility (rock and single-phase fluid) is a constant with a low value In this paper, a novel approach is proposed to perform production data analysis of wells from fractured and fractured-vuggy carbonate reservoirs during the transient flow regime. The proposed rate decline analysis method is very similar to the well test analysis, and the dimensionless variables are based on the conventional well test definitions, as opposed to the Blasingame dimensionless variables. By such, an approach to the rate decline analysis is to reduce the multisolution to certain content. After some reasonable assumption and simplicity, the conventional mathematical modeling of a well in a dual or triple porosity reservoir has been established in a radial cylindrical system [13,16,20], as shown in Appendix A. To solve the modeling, we use the dimensionless effective radius to establish the differential equations. The dimensionless definitions of all kinds of parameters are shown in Appendix B.
The dimensionless mathematical modeling is as follows: For a dual-porosity reservoir [13,20], For a triple-porosity reservoir [19], Definite conditions are shown by taking the tripleporosity reservoir as an example: Initial condition: The well is assumed to produce at constant rate production, and the inner boundary can be written as The outer boundary is assumed to be closed and given by where p mD , p fD , and p vD are dimensionless pressure of matrix, fracture, and vug systems, respectively; ω m , ω f , and ω v are fluid capacitance coefficient of matrix, fracture, and vug systems, respectively; λ mf is the interporosity flow factor of matrix system to fracture system; λ vf is the interporosity flow factor of vug system to fracture system; r D is the dimensionless radial distance; r eD is the dimensionless radial distance of side external boundary; t D is the dimensionless production time; S is the skin factor. Curve 3. Solution to the Mathematical Model. By the use of Laplace transform and inverse Laplace transform methods, the model can be solved. According to the process of solving the model (see Appendix C), the dimensionless pressure in Laplace space for constant rate production can be expressed by [13,16,20,[40][41][42] For a homogenous reservoir, For a dual-porosity reservoir, For a triple-porosity reservoir, Based on the Duhamel principle [43], the solution in Laplace space for modeling well productivity at constant wellbore pressure production can be obtained by where p fD and q D are dimension pressure and production in Laplace space, respectively; u is the Laplace transform variable; I 0 ðÞ is the modified Bessel function of the first kind, zero order; I 1 ðÞ is the modified Bessel function of the first kind, first order; K 0 ðÞ is the modified Bessel function of the second kind, zero order; K 1 () is the modified Bessel function of the second kind, first order.

New Type Curves
3.1. Description of the New Type Curve Method. In real space, the dimensionless production rate (q D ) can be obtained using Stehfest numerical inversion [44] to convert q D back to q D . Therefore, we can obtain the standard log-log type curves of q D , as shown in Figure 3. It can be seen from Figure 3 that the type curves of different dimensionless radial distance of external boundary (r eD ) overlap during the transient flow regime. In addition, the type curves separated gradually with the increase of r eD . This new type curve matching method not only is accurate and practicable but also reduces the multisolution to a certain extent. Particularly, the production well exhibits a long transient flow regime and a boundarydominated regime could not be observed. In order to develop the new type curves, a function related to the dimensionless production rate is introduced [45].
Equation (14) can be rewritten as According to Equation (15), the standard log-log type curves of 1/f ðq D Þ can be obtained. Figure 4 shows the type curves of 1/f ðq D Þ vs. t D for a homogenous reservoir.
Furthermore, we need to introduce the derivative of the function f ′ ðq D Þ to form the rate decline analysis type curves.
The differential form of Equation (14) is written as The derivative of the function f ′ðq D Þ can be written as Taking Equation (16) into Equation (17), one can obtain where q D is the dimensionless production rate and L −1 ðÞ is the inverse Laplace transformation. According to Equation (18), the standard log-log type curves of 1/f ′ ðq D Þ can be obtained. Figure 5 shows the type curves of 1/f ′ ðq D Þ vs. t D for a homogenous reservoir.
Combining Equation (13) with Equations (15) and (18), we can obtain the standard log-log type curves of rate decline analysis of q D , 1/f ðq D Þ, and 1/f ′ðq D Þ vs. t D . Figure 6 shows the type curves of rate decline for a homogenous reservoir. When the theoretical curve is fitted with the measured curve, the permeability, skin coefficient, and other parameters can be interpreted.

Type Curves and Flow Regime Recognition for Multiple
Media. Based on the solutions of the proposed models, the corresponding type curves for multiple media can be obtained. In this section, type curves for production dynamics are plotted and the flow regimes are analyzed as well. Figure 7 shows the standard type curves of rate decline analysis for a well producing located at a dual-porosity reservoir. An entire transient flow process is clearly shown, and the following four main flow periods can be recognized in Figure 7.
Period I is the skin effect and early transition stage. The curves of q D , 1/f ðq D Þ, and 1/f ′ ðq D Þ vs. t D exhibit a downward line.

Geofluids
Period II is the interporosity flow stage of the matrix system to the fracture system. The curve of 1/f ′ðq D Þ vs. t D assumes a ᴧ-shaped hump shape, which represents the interporosity flow of the matrix system to the fracture system. This period is controlled by the storativity ratio of fracture (ω f ) and the interporosity flow coefficient of the matrix system to the fracture system (λ mf ).
Period III is the entire radial flow stage of the matrix and fracture systems. The curves of 1/f ðq D Þ and 1/f ′ ðq D Þ vs. t D converge to a level straight line.
Period IV is the external boundary response stage. The curves of 1/f ðq D Þ and 1/f ′ ðq D Þ vs. t D ultimately converge to a straight line with a slope of "-1," which is the reflection of the pseudosteady state flow. Figure 8 shows the standard type curves of rate decline analysis for a well producing in a triple-porosity reservoir. The following six main flow periods can be recognized: Period I is the skin effect and early radial flow stage. Period II is the interporosity flow stage of the vug system to the fracture system. The curve of 1/f ′ ðq D Þ vs. t D assumes a ∧-shaped hump shape, which represents the interporosity flow of the vug system to the fracture system. The shape and location of ∧-shaped "hump" are controlled by the storativity ratio of vug (ω v ) and the interporosity flow coefficient of the vug system to the fracture system (λ vf ).
Period III is the entire radial flow region of the matrix and fracture systems. The curves of 1/f ðq D Þ and 1/f ′ ðq D Þ vs. t D converge to a level straight line.

Geofluids
Period IV is the interporosity flow stage of the matrix system to the fracture system. As stated above, this period is controlled by the storativity ratio of fracture (ω f ) and the interporosity flow coefficient of the matrix system to the fracture system (λ mf ).
Period V is the entire radial flow region of the vug, matrix, and fracture systems. The curves of 1/f ðq D Þ and 1/f ′ ðq D Þ vs. t D converge to a level straight line.
Period VI is the external boundary response stage. The curves of 1/f ðq D Þ and 1/f ′ðq D Þ vs. t D ultimately converge to a straight line with a slope of "-1." When comparing the type curves of a dual-porosity reservoir with those of a triple-porosity reservoir, the difference is that the curves of 1/f ′ ðq D Þ vs. t D is M-shaped for a tripleporosity reservoir, with the left ∧-shaped being the reflection of the interporosity flow of the vug system to the fracture sys-tem and the right ∧-shaped being the reflection of the interporosity flow of the matrix system to the fracture system.

Characteristics of Transient Production Dynamics.
Based on the proposed model, the effects of some critical parameters on production performance are analyzed. Figures 9 and 10 illustrate the effects of the storativity ratio (ω v ) and interporosity flow coefficient of the vug system to the fracture system (λ vf ) on the production dynamics of a well in fractured-vuggy carbonate reservoirs, respectively. The ω v represents the relative capacity of fluid stored in the vug system. The bigger ω v is the response of relatively abundant reserves in the vug system. The ω v mainly affects the shape of the left ∧-shaped "hump," as shown in Figure 8; a bigger ω v leads to a higher ∧-shaped hump. The λ vf mainly affects the location of the left ∧-shaped "hump" (Figure 9).

Geofluids
A smaller λ vf leads to the later the time of interporosity as λ vf depicts the starting time of the interporosity flow of the vug system to the fracture system. Figures 11 and 12 show the effects of the storativity ratio (ω m ) and interporosity flow coefficient of the matrix system to the fracture system (λ mf ) on type curves, respectively. The ω m represents the relative capacity of fluid stored in the matrix system. As shown in Figure 10, the ω m mainly affects the shape of the right ∧-shaped "hump." A bigger ω m leads to a higher ∧-shaped hump, which means comparatively abundant reserves in the matrix system. The λ mf mainly affects the location of the right ∧-shaped "hump" (Figure 11). Because λ mf represents the starting time of the interporosity flow of the matrix system to the fracture system, therefore, the smaller the λ mf is, the later the time of interporosity is.
The influence of skin factor on type curves is also analyzed, as shown in Figure 13. It can be seen that a smaller skin factor leads to a higher location of the dimensionless rate curve in the early period. In addition, the smaller the skin factor, the slower the rate declines. This is because the incremental value of skin factor results in the increasing additional filtration resistance and the skin effect transition period will last longer. Moreover, the larger the skin factor, the slower the pressure wave propagates to the external boundary and the faster the rate declines. Figure 14 shows the influence of dimensionless radial distance of external boundary on type curves. One can find that r eD mainly influences the late flow periods. The smaller the r eD , the faster the rate declined after the pressure wave propagated to the closed boundary. This is because the external boundary response would be delayed when the r eD increased.

Geofluids
Therefore, the production drop should decrease to maintain the constant BHP of the production well.

Model Verification and Field Application
In this section, we validated the proposed models using two field cases from fractured carbonate oil reservoir and fractured-vuggy carbonate oil reservoir.

Validation Case 1: A Fractured Carbonate Oil Well.
This field example is a vertical well completed in a fractured car-bonate oil reservoir in the Tarim Basin of China. Table 2 shows the basic formation and well parameters that are crucial for well test interpretation and rate decline analysis. The well was put into production on October 15, 2011. During the initial production stage, the well had a production rate of 183.3 m 3 /d. From October 2011 to March 2015, the total production rate declined from 183.3 m 3 /d to 83.7 m 3 /d. The approximate accumulated oil production rate was 24:54 × 10 4 m 3 . The formation pressure was 56.72 MPa.
We first plotted pressure and the pressure derivative curves using the original buildup test data. Then, we chose  9 Geofluids the well-known commercial software of well test analysis Saphir (Ecrin v4.12) developed by the French company Kappa, to conduct well test interpretation. Based on the real field data, we plotted the curves of q/Δp p , ðq/Δp p Þ i , and ðq/Δp p Þ id vs. t ca . We also chose the well-known commercial software of rate decline analysis Topaze (Ecrin v4.12) to perform rate decline interpretation with Blasingame decline type curves. The matching curves of well testing interpretation and rate decline interpretation using a dual-porosity model in Saphir and Topaze are shown in Figures 15(a) and 15(b), respectively. It can be seen from Figures 15(a) and 15(b) that the matching effects are all satisfactory. The main interpretation parameters are shown in Table 2.
At the same time, we plotted the curves of q, 1/f ðqÞ, and 1/f ′ðqÞ vs. t and the proposed method was used to perform rate decline interpretation of the same well. The matching curves of rate decline interpretation by using a dualporosity model with a closed boundary for side are shown in Figure 14(c). As shown, a good match between the new type curves and field data is obtained. The main interpreted parameters using the new method are also presented in Table 3.
According to the matching results shown in Table 2, there is little difference in interpretation results between well test interpretation and the proposed method. However, there is a big difference in the matching results between well test interpretation and the Blasingame method, especially the characteristic parameters S and ω f in the early transient flow period. This is because the production well shows a long transient flow regime and a boundary-dominated regime could not be observed. Compared with the Blasingame method, better history matching results are observed and more reasonably interpreted parameters are obtained, which

Validation Case 2: A Fractured-Vuggy Carbonate Oil
Well. To test the practical applicability of the proposed method, field data obtained from a vertical well completed in a fractured-vuggy carbonate oil reservoir in the Tarim Basin of China is also analyzed. Table 4 summarizes the basic formation and fluid parameters. The well was put into production on January 14, 2010. The initial production rate was 85.84 m 3 /d. The total production rate declined from 85.84 m 3 /d to 51.23 m 3 /d from January 2011 to January 2016. The approximate accumulated oil production rate was 12:78 × 10 4 m 3 , and the formation pressure was 56.91 MPa. Based on the real field data, we also chose the well-known commercial software of well test analysis Saphir and rate decline analysis Topaze (Ecrin v4.12) to perform the well test interpretation and rate decline interpretation, respectively. The matching curves of well testing interpretation and rate decline interpretation using a triple-porosity model in Saphir and Topaze are shown in Figures 16(a) and 16(b), respectively. It can be seen from Figures 16(a) and 16(b) that both of the matching effects are satisfactory. Table 4 shows the main interpretation parameters. Then, the proposed method is used to perform rate decline interpretation of the same well. The matching curves of rate decline interpretation using the new method are shown in Figure 16(c). A satisfactory matching result is obtained, and the main interpreted parameters are also shown in Table 5.
As can be seen from Table 5, the final matching results between well test interpretation and the new method are less than 3%. However, there is a big difference in the interpretation results between well test interpretation and the Blasingame method, especially the characteristic parameters S and ω v in the early transient flow period. This is because the production well is still in transient flow and the reservoir boundaries have not been reached; Blasingame decline curves should not be expected to be applicable. Therefore, the proposed method can accurately undertake production data analysis of wells from fractured and fractured-vuggy carbonate reservoirs during the transient flow regime.

Summary and Conclusions
The highlight of this paper is a novel approach for production data analysis of wells from fractured and fracturedvuggy carbonate reservoirs during the transient flow regime. The main conclusions are as follows:  11 Geofluids (iii) The storativity ratio of the vug system (ω v ), interporosity flow coefficient of the vug system to the fracture system (λ vf ), the storativity ratio of matrix (ω m Þ, interporosity flow of the matrix system to the fracture system (λ mf ), skin factor (S), and dimensionless radial distance of external boundary (r eD ) plays the most significant role in production performance, which should be all considered in the models  12 Geofluids (iv) Two field cases from the Tarim Basin are used to show the accuracy and practicability of the new method. Compared with the Blasingame method, better history matching results are observed and more reasonable interpreted parameters are obtained, which reveal the reliability of the new method in analyzing production data during transient flow regimes