Parameters evaluation of fault-karst carbonate reservoirs with vertical beads-on-string structure based on bottom-hole pressure: Case studies in Shunbei Oilfield, Tarim Basin of Northwestern China

Tarim Basin newly discovered the fault-karst carbonate reservoirs, which are formed by the large-scale tectonic fault activities and multiple-stage 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 beads-on-string 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 multi-fracture-region multi-cave-region series connection physical model by simplifying vertical beads-on-string structure. We consider four kinds of mediums in the proposed physical model, including large caves, small vugs, high-angle 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 beads-on-string 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 unit-slope pressure derivative for the existence of large caves, which is different from the inter-porosity flow regimes for the existence of the vugs (slope ≠ 1); (c) the gravity effect could lead to unit-slope 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 2-weeks test duration in real application. Finally, two cases from Shunbei Oilfield are interpreted to illustrate the practicability and feasibility of proposed method.

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 q Flow rate, m 3 /s q Fluid density, kg/m 3 g Gravitational acceleration, 9.8 m/s 2 t Time, s t D Dimensionless time 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 x n Fracture storativity ratio in fracture region n, dimensionless k n Vug inter-porosity coefficient in fracture region n, dimensionless a Shape factor s Laplace variable, dimensionless

Introduction
Recent years, deeply buried strata (>6000 m) in Tarim Basin have attracted considerable attention. During Chinese 13th Five-Year Plan period, PetroChina and Sinopec increased investment in exploration and achieved major exploration discoveries in the ultra-deep stratas of Tarim Basin. So far, several billion-ton-level oilfields in Tarim Basin have been discovered and developed one after another, such as Tahe, Tazhong, Halahatang, YingMaiLi, and Shunbei Oilfields [1][2][3][4][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 large-scale tectonic fault activities and multiple-stage karstification. Their buried depths are generally more than 6 km [6][7][8]. The ground outcrop, seismic reflection, core sample and well logging data show that the large caves, small vugs, high-angle tectonic fracture and matrix mediums coexist in the fault-karst carbonate reservoirs, causing the complicated flow behaviors (Fig. 1). The rock matrix has extremely low porosity and permeability due to extra-large 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][16][17]. The vugs with low permeability and high storage capacity serve as storage space [18][19][20]. Additionally, the connection patterns of natural fractures and large caves are extremely complex. The vertical beads-on-string structure is regarded as one of the most common patterns in fault-karst 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 beads-on-string structures. Different beads-on-string structures are separated by tight rock matrix. Therefore, each beads-on-string 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 multi-porosity continuum theories are not applicable to accurately describe the complicated flow behaviors [22].
However, the multi-porosity continuum models and radial composite model are often applied to well test interpretation in the early stage of oilfield development [23][24][25][26][27][28][29][30][31][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 dual-porosity 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 non-Darcy 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 fracture-vug horizontal combination in fractured vuggy carbonate reservoirs. Li et al. [10] derived a semianalytical solution for horizontal beads-on-string 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 multi-zonal composite model to describe the flow in horizontal beads-on-string 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 beads-on-string 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 fault-karst carbonate reservoirs with vertical beads-on-string 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 beads-on-string structure. We hope to add a well test analysis method to interpret the measured pressure data in fault-karst carbonate reservoirs. To accomplish this goal, the vertical beadson-string structure is first equivalent to a physical model of large caves vertically interconnecting with high-angle tectonic fractures in series. We consider four types of mediums in the proposed physical model, including large caves, small vugs, high-angle 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.

Physical model
The high-angle tectonic fractures interconnect with large caves in series to form the vertical beads-on-string structure, which is the most common connection pattern in fault-karst carbonate reservoirs. In this subsection, the physical model of cave regions interconnecting with fracture regions in series is presented to characterize the vertical beads-on-string 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 high-angle fractures, vugs and rock matrix. Considering that the permeability of rock matrix is much smaller and the high-angle fracture develops, the main flow in fracture region n is treated as 1-dimensional linear flow along the high-angle fractures and the fluids in the vugs transfer to fracture by pseudo-steady-state 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.
Other assumptions for facilitate the model establishment are given as follows: 1. The vertical beads-on-string 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. 2. 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 . 3. In fracture region n, the high-angle fracture system and the vug system have permeability k fn , k vn , porosity u fn , u vn , comprehensive compressibility C tfn , C tvn and the viscosity l fn , l vn , respectively. The Warren-Root model is used to describe the pseudo-steadystate inter-porosity flow from vugs to fractures. 4. The gravity effect is taken into account due to the flow along z axis.

Mathematical model
Taking two-region-one-cave structure as example, this subsection gives detailed modeling process. As for the multi-region-multi-cave structure, the model establishment and solving steps are given in Appendix A.

Two-region-one-cave model establishment
According to the basic assumption, the main flow in fracture region is 1-D 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 where q is the fluid density, g is gravitational acceleration. Due to the mass transfer from the vugs to the fractures by pseudo-steady-state flow, there is a source term in the mass conservation equation which is given by where q represents the inter-porosity flow volume per unit of time from vug to fracture and Warren and Root model [37] is used to give the inter-porosity 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 For fracture medium; In fracture region I and fracture region II, the initial condition is The inner boundary condition is 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 and with D n (s) is given by where the mobility ratio of fracture region I and fracture region n is defined as M 1n ¼ k , 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 and the corresponding outer boundary condition becomes as follows:

Analytical solution in Laplace space
Solving equations (13) and (14), the general solutions are obtained firstly as follows: where 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 tworegion-one-cave 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].

Type curves
Based on the proposed multi-region-multi-cave model for fault-karst carbonate reservoir with vertical beads-onstring structure, the type curves of two-region-one-cave structure and three-region-two-cave structure are plotted on Figures 3 and 4, respectively. For the case of tworegion-one-cave structure, five typical flow regimes can be recognized on the log-log curves.
1. The first regime is the wellbore-storage regime. The typical characteristic is that the slope of pressure and its derivative curves equals to 1. 2. 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. 3. The third regime is the storage regime of cave I with the typical characteristic of the unit-slope pressure derivative. According to this unique characteristic, the large cave can be identified. 4. The fourth regime is the inter-porosity flow regime, which represents the inter-porosity flow from the vug to fracture nedium. The typical characteristic is that the pressure derivative curves show "V-shape" 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 inter-porosity regime does not equal to 1. Therefore, the large cave and vug mediums can be distinguished by the slope of pressure derivative curves. 5. 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. Figure 4 depicts the type curves of three-region-twocave 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 inter-porosity 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 multi-cave 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 2-weeks 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 long-time production.

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 inter-porosity coefficient, fracture storativity ratio.

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.

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 "V-shape" characteristic at the cave storage regime. With the larger cave I storage coefficient, the "V-shape" characteristic will be wider and deeper.

The effect of gravity coefficient (G D )
Since the vertical beads-on-string structure develops many high-angle 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 unit-slope. 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.

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 inter-porosity 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.

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 "V-shape" characteristic at the inter-porosity 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 "V-shape" characteristic will be wider and deeper.
3. 6 The effect of vug inter-porosity coefficient (k 1 ) In Figure 11, we plot the effects of dimensionless vug interporosity coefficient on type curves (k 1 = 0.0001, k 1 = 0.00001, k 1 = 0.000001). It can be clearly seen that the vug inter-porosity coefficient affects the appearance time of inter-porosity flow regime. With the smaller k 1 , the   inter-porosity flow regime occurs later. This conclusion is similar to the conventional knowledge in multi-porosity continuum theory.

3.7
The effect of fracture storativity ratio (x 1 ) In Figure 12, we further investigate the effects of dimensionless fracture storativity ratio on type curves (x 1 = 0.01, x 1 = 0.03, x 1 = 0.05). Similar to the conventional knowledge, the fracture storativity ratio affects the depth of inter-porosity flow regime. With the smaller fracture storage ratio, the "V-shape" characteristic at inter-porosity flow regime is deeper.

Case studies in Shunbei Oilfield
The Shunbei Oilfield is a typical fault-karst 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, high-angle tectonic fracture and rock matrix coexist in the Shunbei Oilfield. Seismic profile indicates that the vertical beadson-string 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 beads-on-string structure. Moreover, we can clearly observe that the well test curves of SHB A exhibit the cave storage regime with the unit-slope pressure-derivative feature and the inter-porosity flow regime. Therefore, we use the two-region-one-cave model to interpret measured pressure data. The reservoir and well parameters are listed in Table 1. 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.

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 beads-on-string structure. As shown in Figure 15, the seismic profile proves that the SHB B is drilled in the vertical beads-on-string structure. From Figure 16, the well test curves of SHB B exhibit the cave storage regime with the unit-slope pressure derivative feature and the interporosity flow regime. Therefore, we use the proposed tworegion-one-cave model to interpret the measured pressure data. The basic reservoir and well parameters are listed in Table 3.
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.

Conclusion
In this paper, we have presented the multi-region-multicave model for pressure transient analysis of vertical beads-on-string structure, in which four types of mediums including large caves, small vugs, high-angle tectonic fracture and rock matrix are considered. Based on this work, the following conclusions can be reached: 1. For the fault-karst carbonate reservoirs with vertical beads-on-string structure, the typical feature on well test curves is that the cave storage regimes and linear flow regimes alternately appear. 2. The type curves all show the "V-shape" characteristic for the existence of cave and vug mediums. The difference is that the type curve for the existence of large caves will exhibit unit-slope 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. 3. In multi-cave cases, the inter-porosity 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. 4. 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 long-time 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 2-weeks test duration. 5. 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).  e r 1 z 1D e r 2 z 1D Àe r 3 z 2D À e r 4 z 2D r 1 e r 1 z 1D r 2 e r 2 z 1D F 1 : : F 2 : : : : e r 2nÀ1 z ð2nÀ1ÞD r 2nÀ1 e r 2nÀ1 z ð2nÀ1ÞD : : e r 2n z ð2nÀ1ÞD r 2n e r 2n z ð2nÀ1ÞD Àe r 2nþ1 z 2nD F 2nÀ1