Regular Article
Parameters evaluation of faultkarst carbonate reservoirs with vertical beadsonstring structure based on bottomhole pressure: Case studies in Shunbei Oilfield, Tarim Basin of Northwestern China
^{1}
State Key Laboratory of Petroleum Resources and Engineering, China University of PetroleumBeijing, Beijing 102249, PR China
^{2}
Wuqi Oil Production Plant of Yanchang Oilfield Co. Ltd., Wuqi Shaanxi 717600, PR China
^{3}
China National Aviation Fuel Group Logistic Co. Ltd., Shanghai 201100, PR China
^{*} Corresponding authors: chengsq973@163.com; petroyang@163.com
Received:
25
March
2021
Accepted:
21
June
2021
Tarim Basin newly discovered the faultkarst carbonate reservoirs, which are formed by the largescale tectonic fault activities and multiplestage karstification. Four kinds of mediums coexist in the reservoirs, including the large cave, vug, tectonic fracture and matrix. The tectonic fractures interconnect with large caves in series to form the vertical beadsonstring structure, which is the most common connection pattern in reservoirs. To provide a well test method for evaluating this type of structure, this work firstly presents a multifractureregion multicaveregion series connection physical model by simplifying vertical beadsonstring structure. We consider four kinds of mediums in the proposed physical model, including large caves, small vugs, highangle tectonic fracture and rock matrix. The fracture regions mainly contain fracture, vug and matrix mediums. The cave regions contain cave medium. The corresponding mathematical model is also developed, in which the flow in fracture regions obeys the Darcy’s law, while the flow in cave regions is assumed to obey free flow. Furthermore, the gravity is taken into account because the flow is along the vertical direction. Then the typical flow regimes are analyzed and sensitivity analysis is conducted on crucial parameters. Results indicate that (a) the typical feature of vertical beadsonstring structure on type curves is that the cave storage regimes and linear flow regimes alternately appear; (b) the type curves will exhibit the cave storage regimes with unitslope pressure derivative for the existence of large caves, which is different from the interporosity flow regimes for the existence of the vugs (slope ≠ 1); (c) the gravity effect could lead to unitslope pressure and pressure derivative curves, which can be regarded as closed boundary in a peculiar sense; (d) gravity effect is difficult to be observed from well test curves with about 2weeks test duration in real application. Finally, two cases from Shunbei Oilfield are interpreted to illustrate the practicability and feasibility of proposed method.
© C. Wei et al., published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
r_{f}: Equivalent radius of caves and fracture regions, m
C: Wellbore storage coefficient, m^{3}/Pa
C_{vn}: Cave n storage coefficient, m^{3}/Pa
c_{t}: Total system compressibility, Pa^{−1}
c_{f}: Fluid compressibility, Pa^{−1}
c_{φ}: Rock compressibility, Pa^{−1}
h_{2n−1} : Lower boundary of fracture region n, m
h_{2(n−1)} : Upper boundary of fracture region n, m
h_{2n−1} : Upper boundary of cave region n, m
h_{2n} : Lower boundary of cave region n, m
P_{fn}: Fracture pressure in fracture region n, Pa
P_{vn}: Vug pressure in fracture region n, Pa
P_{cvn}: Cave pressure in cave region n, Pa
P_{i}: Initial reservoir pressure, Pa
g: Gravitational acceleration, 9.8 m/s^{2}
G_{D}: Dimensionless gravity coefficient
z_{D}: Dimensionless distance along z axis
C_{D}: Dimensionless wellbore storage coefficient
C_{vnD}: Dimensionless cave n storage coefficient
P_{wfD,CS}: Dimensionless bottomhole pressure with wellbore storage and skin
P_{wfD}: Dimensionless bottomhole pressure
P_{fnD}: Dimensionless fracture pressure in fracture region n
P_{vnD}: Dimensionless vug pressure in fracture region n
P_{cvnD}: Dimensionless cave pressure in cave region n
W_{1n} : The storage ratio of region I to region n, dimensionless
M_{1n} : The mobility ratio of region I to region n, dimensionless
ω_{n}: Fracture storativity ratio in fracture region n, dimensionless
λ_{n}: Vug interporosity coefficient in fracture region n, dimensionless
s: Laplace variable, dimensionless
Subscripts
1 Introduction
Recent years, deeply buried strata (>6000 m) in Tarim Basin have attracted considerable attention. During Chinese 13th FiveYear Plan period, PetroChina and Sinopec increased investment in exploration and achieved major exploration discoveries in the ultradeep stratas of Tarim Basin. So far, several billiontonlevel oilfields in Tarim Basin have been discovered and developed one after another, such as Tahe, Tazhong, Halahatang, YingMaiLi, and Shunbei Oilfields [1–5].
Taking Shunbei Oilfield as example, this type of newly discovered reservoirs behave some unique features that are different from conventional carbonate reservoirs. The reservoirs are formed by largescale tectonic fault activities and multiplestage karstification. Their buried depths are generally more than 6 km [6–8]. The ground outcrop, seismic reflection, core sample and well logging data show that the large caves, small vugs, highangle tectonic fracture and matrix mediums coexist in the faultkarst carbonate reservoirs, causing the complicated flow behaviors (Fig. 1). The rock matrix has extremely low porosity and permeability due to extralarge buried depths, and hence the flow is neglected in most cases [9, 10]. The large caves are usually not filled, whose size can range from meters to even tens of meters [11, 12]. Thus, large caves not only have high storage capacity, but high permeability. The flow inside the large caves no longer obeys the Darcy’s law, but is subjected to free flow [13, 14]. The tectonic fractures with high permeability and low storage capacity serve as main flow medium [15–17]. The vugs with low permeability and high storage capacity serve as storage space [18–20]. Additionally, the connection patterns of natural fractures and large caves are extremely complex. The vertical beadsonstring structure is regarded as one of the most common patterns in faultkarst carbonate reservoirs [5, 21]. As shown in Figure 1a, large caves are widely interconnected with many tectonic fractures in series, which resembles a string of beads and is thus defined as the beadsonstring structures. Different beadsonstring structures are separated by tight rock matrix. Therefore, each beadsonstring structure has its own independent pressure and temperature system, which can be regarded as a “tiny reservoir” in a peculiar sense. Due to these uniqueness, it’s obvious that conventional multiporosity continuum theories are not applicable to accurately describe the complicated flow behaviors [22].
Fig. 1 Geological situations of Shunbei Oilfield (a) vertical beadsonstring structure (Lu et al. [5]); (b) ground outcrop surveying (large cave); (c) rock core sample (vug); (d) well logging results. 
However, the multiporosity continuum models and radial composite model are often applied to well test interpretation in the early stage of oilfield development [23–32]. Oilfield engineers quickly found huge derivation between the interpretation results and actual oilfield conditions. To order to tackle the shortcomings of traditional models, several attempts have been made to evaluate reservoir parameters and well performance. Yang et al. [33] developed a composite model for the case of wells drilled into large caves where the cave is regarded as an equipotential body and the outer fractured region is represented by the dualporosity continuum model. Gao et al. [34] established a similar model where the caves are considered to be filled with loose material and the Barree–Conway model is used to describe the nonDarcy flow in the cave. Wan et al. [14] developed a numerical well test model for wells drilled into filled/unfilled caves and interpreted the measured data in Tarim Basin. Although their models all take into account the existence of large karst caves, they ignore the connection patterns between karst caves and tectonic fractures. Afterwards, Liu et al. [35] established a well test model of multi fracturevug horizontal combination in fractured vuggy carbonate reservoirs. Li et al. [10] derived a semianalytical solution for horizontal beadsonstring fracturecaved carbonate reservoirs. The dimensionless type curves for rate transient analysis are plotted and a case study is conducted to evaluate the reservoir and well properties. Du et al. [12] presented an analytical multizonal composite model to describe the flow in horizontal beadsonstring carbonate reservoirs, which provided a method of determining the cave volume by terms of transient pressure data. Recently, Li et al. [36] developed a well test model for vertical beadsonstring structure. However, they neglected the existence of the vug medium in developing mathematical model. Overall, some oilfield cases cannot obtain satisfactory interpretation results using above analysis methods. Specially, it still lacks pressure analysis approach to characterize the faultkarst carbonate reservoirs with vertical beadsonstring structure.
As a complement to existing models [10, 12, 35, 36], this paper aims at developing a novel well test model for faultkarst carbonate reservoirs with vertical beadsonstring structure. We hope to add a well test analysis method to interpret the measured pressure data in faultkarst carbonate reservoirs. To accomplish this goal, the vertical beadsonstring structure is first equivalent to a physical model of large caves vertically interconnecting with highangle tectonic fractures in series. We consider four types of mediums in the proposed physical model, including large caves, small vugs, highangle tectonic fracture and rock matrix. Following that, corresponding mathematical model is established. Our model also takes the gravity effect into account because the flow is along vertical direction. Then the type curves are analyzed in detail and the sensitivity analysis is conducted on different parameters. Finally, two cases in Shunbei Oilfield are interpreted, which are used to verify the practicability and feasibility of proposed methods.
2 Methodology
2.1 Physical model
The highangle tectonic fractures interconnect with large caves in series to form the vertical beadsonstring structure, which is the most common connection pattern in faultkarst carbonate reservoirs. In this subsection, the physical model of cave regions interconnecting with fracture regions in series is presented to characterize the vertical beadsonstring structure. As shown in Figure 2, the proposed physical model consists of n + 1 fracture regions and n cave regions. The fracture regions and cave regions are all equivalent to cylinders with the radius r_{f}. The wellbore is directly connected with fracture region I. The distance of upper boundary and lower boundary of fracture region n is h_{2(n−1)} and h_{2n−1}, respectively. The fracture region mainly contains highangle fractures, vugs and rock matrix. Considering that the permeability of rock matrix is much smaller and the highangle fracture develops, the main flow in fracture region n is treated as 1dimensional linear flow along the highangle fractures and the fluids in the vugs transfer to fracture by pseudosteadystate interporosity flow. Additionally, the cave region n only contains the large cave, whose upper boundary and lower boundary are h_{2n−1} and h_{2n}, respectively. Two cave regions are connected by a fracture region.
Fig. 2 Schematic of vertical beadsonstring structure in faultkarst carbonate reservoirs (a) ground outcrop surveying; (b) seismic profile reflection; (c) corresponding physical model proposed in this paper. 
Other assumptions for facilitate the model establishment are given as follows:
The vertical beadsonstring structure is filled with single phase and slightly compressible fluid. The initial pressure is P_{i}. The oil well produces at a constant flowrate q.
As the previous literature have mentioned, the pressure drop when the fluid flows across the large caves is very small, which can be neglected [12, 14]. Therefore, each cave region is regarded as a equipotential cylinder, which obeys the free flow. The storage coefficient of large cave n is C_{vn}.
In fracture region n, the highangle fracture system and the vug system have permeability k_{fn}, k_{vn}, porosity φ_{fn}, φ_{vn}, comprehensive compressibility C_{tfn}, C_{tvn} and the viscosity μ_{fn}, μ_{vn}, respectively. The Warren–Root model is used to describe the pseudosteadystate interporosity flow from vugs to fractures.
The gravity effect is taken into account due to the flow along z axis.
2.2 Mathematical model
Taking tworegiononecave structure as example, this subsection gives detailed modeling process. As for the multiregionmulticave structure, the model establishment and solving steps are given in Appendix A.
2.2.1 Tworegiononecave model establishment
According to the basic assumption, the main flow in fracture region is 1D linear flow along the z axis. Since the direction of gravity is opposite with flow velocity, the flow velocity of arbitrary point in the fracture region I is expressed as
(1)where ρ is the fluid density, g is gravitational acceleration.
Due to the mass transfer from the vugs to the fractures by pseudosteadystate flow, there is a source term in the mass conservation equation which is given by
(2)where q represents the interporosity flow volume per unit of time from vug to fracture and Warren and Root model [37] is used to give the interporosity flow flux:
According to equations (1)–(3), the flow governing equation for fracture medium in fracture region I can be written as
The flow governing equation for vug medium in fracture region I is
In the same manner, the flow governing equations in fracture region II are easily given by
In fracture region I and fracture region II, the initial condition is
The inner boundary condition is
(8)and the outer boundary condition is
As for the cave region I, the flow in large caves obeys the free flow. It means that the pressure spreads fastly across the large cave. According to the previous work [12], the pressure drop in the large cave can be neglected. Therefore, the cave region I is regarded as a equipotential body and the pressure continuity condition at the interface can be given by
Additionally, there are three reasons that will cause the fluid mass in the cave region I to change.
The oil in the cave region I flows to the fracture region I.
The oil transfers from the fracture region II to the cave region I.
The fluid mass changes due to the compressibility of fluid and large cave, which is characterized by the cave I storage.
According to the mass conservation law, the flowrate continuity condition at the cave region I can be described as
The initial condition in the cave region I is
To facilitate solving the established model, we need to transform all equations dimensionlessly. The dimensionless definitions are given in Table B1 (see Appendix B). Applying the dimensionless and Laplace transformation to equations (4)–(12), we arrive to the following equations.
In fracture regions I and II, the flow governing equations are transformed into
(15)where the mobility ratio of fracture region I and fracture region n is defined as , the storage ratio of fracture region I and fracture region n is defined as .
In cave region I, the pressure and flowrate continuity conditions are transformed into
The initial conditions of the fracture regions I, II and cave region I become
The inner boundary condition becomes
(19)and the corresponding outer boundary condition becomes as follows:
2.2.2 Analytical solution in Laplace space
Solving equations (13) and (14), the general solutions are obtained firstly as follows:
Then the c_{1}–c_{4} is solved by means of equations (16)–(20):
At z = 0, the bottomhole pressure equals to P_{f1D}. According to equation (21), the bottomhole pressure solution of tworegiononecave model is given by
If the wellbore storage and skin are considered, the Van Everdingen and Hurst [38] work can be refered:
Finally, P_{wfD,CS} is obtained by the Stehfest numerical inversion [39].
2.3 Type curves
Based on the proposed multiregionmulticave model for faultkarst carbonate reservoir with vertical beadsonstring structure, the type curves of tworegiononecave structure and threeregiontwocave structure are plotted on Figures 3 and 4, respectively. For the case of tworegiononecave structure, five typical flow regimes can be recognized on the loglog curves.
The first regime is the wellborestorage regime. The typical characteristic is that the slope of pressure and its derivative curves equals to 1.
The second regime is the linear flow regime, which represents the flow in fracture region I. The typical characteristic is that the slope of pressure and its derivative curves equals to 1/2.
The third regime is the storage regime of cave I with the typical characteristic of the unitslope pressure derivative. According to this unique characteristic, the large cave can be identified.
The fourth regime is the interporosity flow regime, which represents the interporosity flow from the vug to fracture nedium. The typical characteristic is that the pressure derivative curves show “Vshape” characteristic, which is similar with the characteristic of large cave. In fact, the typical characteristics of large cave and vug are different on well test curves. As shown in Figure 3, the slope of pressure derivative at interporosity regime does not equal to 1. Therefore, the large cave and vug mediums can be distinguished by the slope of pressure derivative curves.
The fifth regime is boundary control flow regime. The typical characteristic is that the slope of pressure and its derivative curves equals to 1/2, which represents the linear flow in fracture region. However, this characteristic is different from that of conventional infinite boundary. In conventional models with infinite boundary, the pressure derivative curves are 0.5 horizontal lines.
Fig. 3 Type curves of tworegiononecave structure. 
Fig. 4 Type curves of threeregiontwocave structure. 
Figure 4 depicts the type curves of threeregiontwocave structure. In this case, the demensionless wellbore storage and skin are all set to 0. The type curves contains three linear flow regimes and two cave storage regimes, corresponding to three fracture regions and two cave regions. Therefore, the number of cave storage regime is corresponding to the number of large caves. Additionally, no interporosity flow regime is observed on type curves because these regime are masked by the cave storage regimes, causing that the vug medium usually cannot be identified by the pressure and its derivative curves in multicave cases. At this point, other information (e.g., the core sample, well logging and ground outcrop surveying) should be combined to more accurately characterize the reservoirs.
This paper also takes the gravity effect into account. As shown in Figure 4, the gravity effect results in similar flow behaviors with the closed boundary that the pressure and derivative curves go upward with unit slope at late times. We further evaluate the appearance time of gravity effect. According to the definition of dimensionless gravity coefficient, the dimensionless gravity coefficient of crude oil is between 10^{−3} and 10^{−4}. We set three different G_{D} to analyze the appearance time of gravity effect (G_{D} = 0.001, G_{D} = 0.0005, G_{D} = 0.0001). As can be clearly seen from Figure 5, the appearance time of gravity effect is late (about t_{D} = 3 × 10^{6}). It’s a difficult thing that the gravity effect can be observed from well test curves with about 2weeks test duration in real application. However, this phenomenon gives the engineers a consideration that necessary treatments (e.g., water injection) should be taken to offset the gravity effects during longtime production.
Fig. 5 Appearance time of gravity effect. 
3 Sensitivity analysis
To get a better understanding of the proposed model, this subsection analyzes the effect of several parameters on the type curves, including the fracture region I length, cave I storage coefficient, gravity coefficient, mobility ratio of region I to region II, storage ratio of region I to region II, vug interporosity coefficient, fracture storativity ratio.
3.1 The effect of fracture region I length (z_{1D})
To investigate the effect of dimensionless fracture region I length (Z_{1D}) on type curves, we conduct three analysis cases (Z_{1D} = 3, Z_{1D} = 4, Z_{1D} = 5). As shown in Figure 6, the fracture region I length affects the duration time of linear flow regime and the appearance time of cave storage regime. With the larger Z_{1D}, the linear flow regime lasts longer and the cave storage regime appears later. The reason for this is that the larger Z_{1D} means the longer fracture region. Pressure wave needs to take more time to pass through the fracture region.
Fig. 6 The effect of dimensionless fracture region I length (z_{1D}) on type curves. 
3.2 The effect of cave I storage coefficient (C_{v1D})
Three cases of different dimensionless cave I storage coefficient (C_{v1D} = 5, C_{v1D} = 10, C_{v1D} = 15) are designed to investigate their effects on type curves. As shown in Figure 7, the C_{v1D} affects the depth and width of “Vshape” characteristic at the cave storage regime. With the larger cave I storage coefficient, the “Vshape” characteristic will be wider and deeper.
Fig. 7 The effect of dimensionless cave I storage coefficient (C_{v1D}) on type curves. 
3.3 The effect of gravity coefficient (G_{D})
Since the vertical beadsonstring structure develops many highangle tectonic fractures and the main flow is along vertical direction, the gravity effect should be a very important factor, deserving to be discussed. This part designs three cases (G_{D} = 0, G_{D} = 0.002, G_{D} = 0.005) to investigate the effect of dimensionless gravity coefficient. As shown in Figure 8, the gravity effect behaves similar with that of closed boundary. At late time, the pressure and its derivative go upward with unitslope. With the larger dimensionless gravity coefficient, the pressure and its derivative curves go upward earlier. Therefore, the gravity factor can be seemed as closed boundary in a peculiar sense. The gravity effect appears earlier with the larger gravity. At this point, the magnitude of gravity can be equivalent to the distance of closed boundary.
Fig. 8 The effect of dimensionless gravity coefficient (G_{D}) on type curves. 
3.4 The effect of mobility ratio of region I to region II (M_{12})
Three cases of different dimensionless mobility ratio of region I to region II (M_{12} = 0.5, M_{12} = 1, M_{12} = 2) are designed to investigate their effects on type curves. As shown in Figure 9, it can be clearly observed that the larger the M_{12} is, the later the time of interporosity flow regime appears. The reason for this is that the larger M_{12} means the lower mobility in fracture region II, leading to the later appearance time of flow regimes in fracture region II.
Fig. 9 The effect of dimensionless mobility ratio of region I to region II (M_{12}) on type curves. 
3.5 The effect of storage ratio of region I to region II (W_{12})
To further investigate the effect of dimensionless storage ratio of region I to region II (W_{12}), three cases (W_{12} = 2, W_{12} = 1, W_{12} = 0.5) are designed to analyze their effects on type curves. It can be clearly seen from Figure 10 that the “Vshape” characteristic at the interporosity flow regime will be wider and deeper with the smaller W_{12}. It’s because the smaller W_{12} represents the bigger storage capacity of fracture region II. Therefore, the “Vshape” characteristic will be wider and deeper.
Fig. 10 The effect of dimensionless storage ratio of region I to region II (W_{12}) on type curves. 
3.6 The effect of vug interporosity coefficient (λ_{1})
In Figure 11, we plot the effects of dimensionless vug interporosity coefficient on type curves (λ_{1} = 0.0001, λ_{1} = 0.00001, λ_{1} = 0.000001). It can be clearly seen that the vug interporosity coefficient affects the appearance time of interporosity flow regime. With the smaller λ_{1}, the interporosity flow regime occurs later. This conclusion is similar to the conventional knowledge in multiporosity continuum theory.
Fig. 11 The effect of dimensionless vug interporosity coefficient (λ_{1}) on type curves. 
3.7 The effect of fracture storativity ratio (ω_{1})
In Figure 12, we further investigate the effects of dimensionless fracture storativity ratio on type curves (ω_{1} = 0.01, ω_{1} = 0.03, ω_{1} = 0.05). Similar to the conventional knowledge, the fracture storativity ratio affects the depth of interporosity flow regime. With the smaller fracture storage ratio, the “Vshape” characteristic at interporosity flow regime is deeper.
Fig. 12 The effect of dimensionless fracture storativity ratio (ω_{1}) on type curves. 
4 Case studies in Shunbei Oilfield
The Shunbei Oilfield is a typical faultkarst carbonate reservoir in Tarim Basin of Northwestern China, which is put into production in recent years. The ground outcrop surveying, drilling information, core sample and well logging data show that large caves, small vugs, highangle tectonic fracture and rock matrix coexist in the Shunbei Oilfield. Seismic profile indicates that the vertical beadsonstring structure is one of the most common connectivity patterns of tectonic fractures and large caves. To illustrate the practicability and feasibility, we apply the proposed model to interpret two cases chosen from Shunbei Oilfield.
Field case 1
The well SHB A produces at a constant production (65 m^{3}/d). We choose the SHB A as a case study, which is conducted a buildup test for 255 h. As shown in Figure 13, the seismic profile indicates that the SHB A is drilled in the vertical beadsonstring structure. Moreover, we can clearly observe that the well test curves of SHB A exhibit the cave storage regime with the unitslope pressurederivative feature and the interporosity flow regime. Therefore, we use the tworegiononecave model to interpret measured pressure data. The reservoir and well parameters are listed in Table 1.
Fig. 13 Seismic profile reflection of Well SHB A. 
Basic parameters of Well SHB A.
As shown in Figure 14, the results of theoretical model achieve good fitting quality with the recorded pressure and derivative data. The interpretation results are given in Table 2. Besides the conventional parameters, the equivalent cave radius and cave height are obtained by using the proposed model. Based on these results, the equivalent cave volume can be further calculated.
Fig. 14 The application of the proposed model on Well SHB A. 
The interpretations results of Well SHB A.
Field case 2
The well SHB B produces at a constant production (60 m^{3}/d), which is conducted a buildup test for 231 h. The seismic profile proves that the SHB B is drilled in the vertical beadsonstring structure. As shown in Figure 15, the seismic profile proves that the SHB B is drilled in the vertical beadsonstring structure. From Figure 16, the well test curves of SHB B exhibit the cave storage regime with the unitslope pressure derivative feature and the interporosity flow regime.. Therefore, we use the proposed tworegiononecave model to interpret the measured pressure data. The basic reservoir and well parameters are listed in Table 3.
Fig. 15 Seismic profile reflection of Well SHB B. 
Fig. 16 The application of the proposed model on Well SHB B. 
Basic parameters of Well SHB B.
As shown in Figure 16, the results of theoretical model show good fitting with the recorded pressure and derivative data. The interpretation results are given in Table 4.
The interpretation results of Well SHB B.
5 Conclusion
In this paper, we have presented the multiregionmulticave model for pressure transient analysis of vertical beadsonstring structure, in which four types of mediums including large caves, small vugs, highangle tectonic fracture and rock matrix are considered. Based on this work, the following conclusions can be reached:
For the faultkarst carbonate reservoirs with vertical beadsonstring structure, the typical feature on well test curves is that the cave storage regimes and linear flow regimes alternately appear.
The type curves all show the “Vshape” characteristic for the existence of cave and vug mediums. The difference is that the type curve for the existence of large caves will exhibit unitslope pressure derivative, while this slope for vugs is not equal to 1. Therefore, cave and vug medium can be distinguished by the slope of pressure derivative curves.
In multicave cases, the interporosity flow regime usually is masked by the cave storage regimes, causing that the vugs cannot be identified from type curves. At this point, other information (e.g., core sample, well logging) should be combined to more accurately characterize reservoirs.
The gravity effect could result in similar pressure behaviors with closed boundary that the pressure and its derivative curves go upward with unit slope at late times. This phenomenon gives the engineers a consideration that necessary treatments should be taken to offset the gravity effects during longtime production. Moreover, the gravity effect appear earlier with the larger gravity and it’s difficult to observe the gravity effect from well test curves with about 2weeks test duration.
Two cases from the Shunbei Oilfield are interpreted using the proposed model. Results indicate that our models are feasible for well test analysis. Some key reservoir parameters that are not readily accessible from conventional methods can be obtained (e.g., equivalent cave volume).
Acknowledgments
The authors would like to express our thanks for the financial support from The National Natural Science Fund of China (No. 11872073), the Strategic Cooperation Technology Projects of CNPC and CUPB (ZLZX202002), and the Fundamental Research Funds for the Central Universities (2462020XKBH003, 2462020YXZZ028).
Appendix A – Multiregionmulticave model establishment
The multiregionmulticave model contains n + 1 fracture regions and n cave regions. This section gives three parts to introduce the detailed modeling derivations. The whole process is similar with tworegiononecave model.
A.1 Flow equations and boundary conditions
First, the flow equations of n + 1 fracture regions are given by terms of Warren–Root model:
(A1)and the pressure and flowrate continuity conditions of n cave regions are given as follows:
Then the initial conditions of n + 1 fracture regions and n cave regions are respectively
The inner boundary condition is
(A5)and the outer boundary condition is
A.2 Dimensionless forms in Laplace space
Applying dimensionless transformation and Laplace transformation to equations (A1)–(A6), we arrive to the following equations.
The flow governing equations of n + 1 fracture regions are transformed into
(A7)and the pressure and flowrate continuity conditions of n cave regions are transformed into
The initial conditions of n + 1 fracture regions and n cave regions are transformed into
The inner boundary condition is transformed into
(A11)and the corresponding outer boundary condition is transformed into
A.3 Analytical solution in Laplace space
Solving equation (A7), the general solutions are first obtained:
Equation (A13) contains 2n + 2 unknowns (c_{1} − c_{2n + 2}). According to equations (A8)–(A12), c_{1} − c_{2n + 2} can be uniquely determined. Since there are too many equations, we construct a matrix equation to solve c_{1} − c_{2n + 2}:
In equation (A15), F_{2n − 1} and F_{2n} are defined as follows:
Solving equation (A15), we can get the c_{1 −} c_{2n + 2}. Then the bottomhole pressure is easily given according to equation (24),
If wellbore storage and skin are considered, the Van Everdingen and Hurst [38] work can be refered:
Eventually, we obtain the solutions of multiregionmulticave model.
Appendix B – Dimensionless definitions
Dimensionless definition.
References
 Zhu G., Zou C., Yang H., Wang K., Zheng D., Zhu Y.F., Wang Y. (2015) Hydrocarbon accumulation mechanisms and industrial exploration depth of largearea fracture–cavity carbonates in the Tarim Basin, western China, J. Petrol. Sci. Eng. 133, 889–907. https://doi.org/10.1016/j.petrol.2015.03.014. [CrossRef] [Google Scholar]
 Li Y., Hou J.G., Li Y.Q. (2016) Features and classified hierarchical modeling of carbonate fracturecavity reservoirs, Pet. Explor. Dev. 43, 4, 655–662. [CrossRef] [Google Scholar]
 Zhu G., Wang H., Weng N. (2016) TSRaltered oil with highabundance thiaadamantanes of a deepburied Cambrian gas condensate reservoir in Tarim Basin, Marine Petrol. Geol. 69, 1–12. [CrossRef] [Google Scholar]
 Xiao Z., Li M., Huang S., Wang T., Zhang B., Fang R., Zhang K., Ni Z., Zhao Q., Wang D. (2016) Source, oil charging history and filling pathways of the Ordovician carbonate reservoir in the Halahatang Oilfield, Tarim Basin, NW China, Marine Petrol. Geol. 73, 59–71. [CrossRef] [Google Scholar]
 Lu X., Wang Y., Tian F., Li X., Yang D., Li T., Lv Y., He X. (2017) New insights into the carbonate karstic fault system and reservoir formation in the Southern Tahe area of the Tarim Basin, Marine and Petroleum Geology 86, 587–605. https://doi.org/10.1016/j.marpetgeo.2017.06.023. [CrossRef] [Google Scholar]
 Zhu G., Weng N., Wang H., Yang H., Zhang S., Su J., Liao F., Zhang B., Ji Y. (2015) Origin of diamondoid and sulphur compounds in the Tazhong Ordovician condensate, Tarim Basin, China: Implications for hydrocarbon exploration in deepburied strata, Marine Petrol. Geol. 62, 14–27. [CrossRef] [Google Scholar]
 Jia C., Pang X. (2015) Research processes and main development directions of deep hydrocarbon geological theories, Acta Petrolei Sinica 36, 12, 1457–1469. [Google Scholar]
 Zhang L., Guo X., Hao F., Zou H., Li P. (2016) Lithologic characteristics and diagenesis of the Upper Triassic Xujiahe formation, Yuanba area, northeastern Sichuan Basin, J. Natural Gas Sci. Eng. 35, 1320–1335. [CrossRef] [Google Scholar]
 Zhang K. (2012) Strategic replacement situation and outlook of China oilgas production area, Pet. Explor. Dev. 395, 547–559. [CrossRef] [Google Scholar]
 Li Y., Yu Q., Jia C., Liu P., Wang Q., Wang D. (2020) Rate transient analysis for coupling Darcy flow and free flow in beadstring fracturecaved carbonate reservoirs, J. Pet. Sci. Eng. 195, 107809. https://doi.org/10.1016/j.petrol.2020.107809. [CrossRef] [Google Scholar]
 Zhao W., Shen A., Qiao Z., Zheng J., Wang X. (2014) Carbonate karst reservoirs of the Tarim Basin, northwest China: Types, features, origins, and implications for hydrocarbon exploration, Interpretation 23, SF65SF90. [CrossRef] [Google Scholar]
 Du X., Li Q., Lu Z., Li P., Xian Y., Xu Y., Li D., Lu D. (2020) Pressure transient analysis for multivug composite fractured vuggy carbonate reservoirs, J. Pet. Sci. Eng. 193, 107389. https://doi.org/10.1016/j.petrol.2020.107389. [CrossRef] [Google Scholar]
 Huang Z.Q., Yao J., Wang Y.Y. (2013) An efficient numerical model for immiscible twophase flow in fractured karst reservoirs, Commun. Comput. Phys. 13, 2, 540–558. [CrossRef] [Google Scholar]
 Wan Y.Z., Liu Y.W., Chen F.F., Wu N.Y., Hu G.W. (2018) Numerical well test model for caved carbonate reservoirs and its application in Tarim Basin, China, J. Pet. Sci. Eng. 161, 611–624. https://doi.org/10.1016/j.petrol.2017.12.013. [CrossRef] [Google Scholar]
 Flemisch B., Berre I., Boon W., Fumagalli A., Schwenck N., Scotti A., Stefansson I., Tatomir A. (2018) Benchmarks for singlephase flow in fractured porous media, Adv. Water Resour. 111, 239–258. https://doi.org/10.1016/j.advwatres.2017.10.036. [CrossRef] [Google Scholar]
 Berre I., Doster F., Keilegavlen E. (2019) Flow in fractured porous media: A review of conceptual models and discretization approaches, Transp. Porous Media 130, 1, 215–236. [CrossRef] [Google Scholar]
 Wang D., Sun J., Li Y., Peng H. (2019) An efficient hybrid model for nonlinear twophase flow in fractured lowpermeability reservoir, Energies 12, 15, 2850. [CrossRef] [Google Scholar]
 Wei C., Cheng S., Wang Y., Shi W., Li J., Zhang J., Yu H. (2021) Practical pressuretransient analysis solutions for a well intercepted by finite conductivity vertical fracture in naturally fractured reservoirs, J. Pet. Sci. Eng. 204, 108768. https://doi.org/10.1016/j.petrol.2021.108768. [CrossRef] [Google Scholar]
 Dejam M., Hassanzadeh H., Chen Z. (2018) Semianalytical solution for pressure transient analysis of a hydraulically fractured vertical well in a bounded dualporosity reservoir, J. hydrol. 565, 289–301. https://doi.org/10.1016/j.jhydrol.2018.08.020. [CrossRef] [Google Scholar]
 Yan X., Huang Z., Yao J., Zhang Z., Liu P., Li Y., Fan D. (2019) Numerical simulation of hydromechanical coupling in fractured vuggy porous media using the equivalent continuum model and embedded discrete fracture model, Adv. Water Resour. 126, 137–154. [Google Scholar]
 Tian F., Jin Q., Lu X., Lei Y., Zhang L., Zheng S., Zhang H., Rong Y., Liu N. (2016) Multilayered Ordovician paleokarst reservoir detection and spatial delineation: A case study in the Tahe Oilfield, Tarim Basin, Western China, Marine Petrol. Geol. 69, 53–73. [CrossRef] [Google Scholar]
 Popov P., Qin G., Bi L., Efendiev Y., Ewing R.E., Li J. (2009) Multiphysics and multiscale methods for modeling fluid flow through naturally fractured carbonate karst reservoirs, SPE Reserv. Eval. Eng. 12, 02, 218–231. [CrossRef] [Google Scholar]
 Abdassah D., Ershaghi I. (1986) Tripleporosity systems for representing naturally fractured reservoirs, SPE Form. Eval. 1, 02, 113–127. [CrossRef] [Google Scholar]
 Chu W.C., Shank G.D. (1993) A new model for a fractured well in a radial, composite reservoir (includes associated papers 27919, 28665 and 29212), SPE Form. Eval. 8, 03, 225–232. https://doi.org/10.2118/20579PA. [CrossRef] [Google Scholar]
 Kossack C.A., Gurpinar O. (2001) A methodology for simulation of vuggy and fractured reservoirs, in: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. https://doi.org/10.2118/66366MS. [Google Scholar]
 CamachoVelazquez R., VasquezCruz M., CastrejonAivar R., AranaOrtiz V. (2002) Pressure transient and decline curve behaviors in naturally fractured Vuggy carbonate reservoirs, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
 Wu Y.S., EhligEconomides C.A., Qin G., Kang Z., Zhang W., Ajayi B.T., Tao Q. (2007) A triplecontinuum pressuretransient model for a naturally fractured vuggy reservoir, in: Paper presented at the SPE Annual Technical Conference and Exhibition, Anaheim, California, USA, Society of Petroleum Engineers. https://doi.org/10.2118/110044MS. [Google Scholar]
 Nie R.S., Meng Y.F., Guo J.C., Jia Y.L. (2012) Modeling transient flow behavior of a horizontal well in a coal seam, Int. J. Coal Geol. 92, 54–68. https://doi.org/10.1016/j.coal.2011.12.005. [CrossRef] [Google Scholar]
 Zheng D., Yuan B., Moghanloo R.G. (2017) Analytical modeling dynamic drainage volume for transient flow towards multistage fractured wells in composite shale reservoirs, J. Pet. Sci. Eng. 149, 756–764. https://doi.org/10.1016/j.petrol.2016.11.023. [CrossRef] [Google Scholar]
 Abbasi M., Madani M., Sharifi M., Kazemi A. (2018) Fluid flow in fractured reservoirs: Exact analytical solution for transient dual porosity model with variable rock matrix block size, J. Pet. Sci. Eng. 164, 571–583. https://doi.org/10.1016/j.petrol.2018.01.010. [CrossRef] [Google Scholar]
 Wang Y., Cheng S., Zhang K., An X. (2020) Investigation on the pressure behavior of injectors influenced by waterfloodinduced fractures: field cases in Huaqing reservoir, Changqing Oilfield, China, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 75, 20. [CrossRef] [Google Scholar]
 Shi W., Yao Y., Cheng S., Li H., Wang M., Cui N., Zhang C., Li H., Tu K., Shi Z. (2021) Investigation on the pressure response behavior of twolayer vertical mixed boundary reservoir: field cases in Western Sichuan XC gas field, China, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 76, 2. [CrossRef] [Google Scholar]
 Yang F., Wang X.H., Liu H. (2011) Well test interpretation model for wells drilled in cavity of fractured vuggy carbonate reservoirs, Chin. J. Hydrodyn 26, 278–283. [Google Scholar]
 Gao B., Huang Z.Q., Yao J., Lv X.R., Wu Y.S. (2016) Pressure transient analysis of a well penetrating a filled cavity in naturally fractured carbonate reservoirs, J. Pet. Sci. Eng. 145, 392–403. https://doi.org/10.1016/j.petrol.2016.05.037. [CrossRef] [Google Scholar]
 Liu J.Y., Liu Z.B., Zou N., Liu X.L., You X.T., Yi S. (2020) A well test model study of multi fracturevug combination for fractured vuggy carbonate reservoirs, in: ICAE International Conference on Applied Energy, 2020, December 1–10, 2020, Bangkok. [Google Scholar]
 Li Q., Du X., Tang Q., Xu Y., Li P., Lu D. (2021) A novel well test model for fractured vuggy carbonate reservoirs with the vertical beadonastring structure, J. Pet. Sci. Eng. 196, 107938. https://doi.org/10.1016/j.petrol.2020.107938. [CrossRef] [Google Scholar]
 Warren J.E., Root P.J. (1963) The behavior of naturally fractured reservoirs, Soc. Pet. Eng. J. 3, 03, 245–255. [Google Scholar]
 Van Everdingen A.F., Hurst W. (1949) The application of the Laplace transformation to flow problems in reservoirs, J. Pet. Technol. 1, 12, 305–324. [CrossRef] [Google Scholar]
 Stehfest H. (1970) Algorithm 368: Numerical inversion of Laplace transforms [D5], Commun. ACM 13, 1, 47–49. [Google Scholar]
All Tables
All Figures
Fig. 1 Geological situations of Shunbei Oilfield (a) vertical beadsonstring structure (Lu et al. [5]); (b) ground outcrop surveying (large cave); (c) rock core sample (vug); (d) well logging results. 

In the text 
Fig. 2 Schematic of vertical beadsonstring structure in faultkarst carbonate reservoirs (a) ground outcrop surveying; (b) seismic profile reflection; (c) corresponding physical model proposed in this paper. 

In the text 
Fig. 3 Type curves of tworegiononecave structure. 

In the text 
Fig. 4 Type curves of threeregiontwocave structure. 

In the text 
Fig. 5 Appearance time of gravity effect. 

In the text 
Fig. 6 The effect of dimensionless fracture region I length (z_{1D}) on type curves. 

In the text 
Fig. 7 The effect of dimensionless cave I storage coefficient (C_{v1D}) on type curves. 

In the text 
Fig. 8 The effect of dimensionless gravity coefficient (G_{D}) on type curves. 

In the text 
Fig. 9 The effect of dimensionless mobility ratio of region I to region II (M_{12}) on type curves. 

In the text 
Fig. 10 The effect of dimensionless storage ratio of region I to region II (W_{12}) on type curves. 

In the text 
Fig. 11 The effect of dimensionless vug interporosity coefficient (λ_{1}) on type curves. 

In the text 
Fig. 12 The effect of dimensionless fracture storativity ratio (ω_{1}) on type curves. 

In the text 
Fig. 13 Seismic profile reflection of Well SHB A. 

In the text 
Fig. 14 The application of the proposed model on Well SHB A. 

In the text 
Fig. 15 Seismic profile reflection of Well SHB B. 

In the text 
Fig. 16 The application of the proposed model on Well SHB B. 

In the text 