© C. Wei et al., published by IFP Energies nouvelles, 2021

Nomenclature

r_f: Equivalent radius of caves and fracture regions, m

z: Distance along z axis, m

φ: Porosity, fraction

μ: Viscosity, Pa s

B: Volume factor

C: Wellbore storage coefficient, m³/Pa

C_vn: Cave n storage coefficient, m³/Pa

S: Skin factor

c_t: Total system compressibility, Pa⁻¹

c_f: Fluid compressibility, Pa⁻¹

c_φ: Rock compressibility, Pa⁻¹

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

q: Flow rate, m³/s

ρ: Fluid density, kg/m³

g: Gravitational acceleration, 9.8 m/s²

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

ω_n: Fracture storativity ratio in fracture region n, dimensionless

λ_n: Vug inter-porosity coefficient in fracture region n, dimensionless

α: Shape factor

s: Laplace variable, dimensionless

Subscripts

i: Initial condition

D: Dimensionless

f: Natural fracture

v: Vugs

cv: Caves

s: Laplace space

1 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–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–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–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 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–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 semi-analytical solution for horizontal beads-on-string fracture-caved 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 fault-karst 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 beads-on-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.

2 Methodology

2.1 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 inter-porosity 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 φ_fn, φ_vn, comprehensive compressibility C_tfn, C_tvn and the viscosity μ_fn, μ_vn, respectively. The Warren–Root model is used to describe the pseudo-steady-state inter-porosity flow from vugs to fractures.

4. The gravity effect is taken into account due to the flow along z axis.

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

2.2.1 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

(1)where ρ 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

(2)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:

(3)

According to equations (1)–(3), the flow governing equation for fracture medium in fracture region I can be written as

(4)

The flow governing equation for vug medium in fracture region I is

(5)

In the same manner, the flow governing equations in fracture region II are easily given by

(6)

In fracture region I and fracture region II, the initial condition is

(7)

The inner boundary condition is

(8)and the outer boundary condition is

(9)

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

(10)

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

(11)

The initial condition in the cave region I is

(12)

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

(13)and

(14)with D_n(s) is given by

(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

(16)

(17)

The initial conditions of the fracture regions I, II and cave region I become

(18)

The inner boundary condition becomes

(19)and the corresponding outer boundary condition becomes as follows:

(20)

2.2.2 Analytical solution in Laplace space

Solving equations (13) and (14), the general solutions are obtained firstly as follows:

(21)

(22)where

Then the c₁–c₄ is solved by means of equations (16)–(20):

(23)

At z = 0, the bottomhole pressure equals to P_f1D. According to equation (21), the bottomhole pressure solution of two-region-one-cave model is given by

(24)

If the wellbore storage and skin are considered, the Van Everdingen and Hurst [38] work can be refered:

(25)

Finally, P_wfD,CS is obtained by the Stehfest numerical inversion [39].

2.3 Type curves

Based on the proposed multi-region-multi-cave model for fault-karst carbonate reservoir with vertical beads-on-string 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 two-region-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-two-cave 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⁻³ and 10⁻⁴. 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⁶). 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.

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

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

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

3.4 The effect of mobility ratio of region I to region II (M₁₂)

Three cases of different dimensionless mobility ratio of region I to region II (M₁₂ = 0.5, M₁₂ = 1, M₁₂ = 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₁₂ is, the later the time of inter-porosity flow regime appears. The reason for this is that the larger M₁₂ means the lower mobility in fracture region II, leading to the later appearance time of flow regimes in fracture region II.

3.5 The effect of storage ratio of region I to region II (W₁₂)

To further investigate the effect of dimensionless storage ratio of region I to region II (W₁₂), three cases (W₁₂ = 2, W₁₂ = 1, W₁₂ = 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₁₂. It’s because the smaller W₁₂ 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 (λ₁)

In Figure 11, we plot the effects of dimensionless vug inter-porosity coefficient on type curves (λ₁ = 0.0001, λ₁ = 0.00001, λ₁ = 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 λ₁, 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 (ω₁)

In Figure 12, we further investigate the effects of dimensionless fracture storativity ratio on type curves (ω₁ = 0.01, ω₁ = 0.03, ω₁ = 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.

4 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 beads-on-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³/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³/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 inter-porosity flow regime.. Therefore, we use the proposed two-region-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.

5 Conclusion

In this paper, we have presented the multi-region-multi-cave 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).

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 (ZLZX2020-02), and the Fundamental Research Funds for the Central Universities (2462020XKBH003, 2462020YXZZ028).

Appendix A – Multi-region-multi-cave model establishment

The multi-region-multi-cave 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 two-region-one-cave 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:

(A2)

Then the initial conditions of n + 1 fracture regions and n cave regions are respectively

(A3)

(A4)

The inner boundary condition is

(A5)and the outer boundary condition is

(A6)

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

(A8)

The initial conditions of n + 1 fracture regions and n cave regions are transformed into

(A9)

(A10)

The inner boundary condition is transformed into

(A11)and the corresponding outer boundary condition is transformed into

(A12)

A.3 Analytical solution in Laplace space

Solving equation (A7), the general solutions are first obtained:

(A13)where

(A14)

Equation (A13) contains 2n + 2 unknowns (c₁ − c_{2n + 2}). According to equations (A8)–(A12), c₁ − c_{2n + 2} can be uniquely determined. Since there are too many equations, we construct a matrix equation to solve c₁ − c_{2n + 2}:

(A15)

In equation (A15), F_{2n − 1} and F_2n are defined as follows:

(A16)

(A17)

Solving equation (A15), we can get the c_{1 −} c_{2n + 2}. Then the bottomhole pressure is easily given according to equation (24),

(A18)

If wellbore storage and skin are considered, the Van Everdingen and Hurst [38] work can be refered:

(A19)

Eventually, we obtain the solutions of multi-region-multi-cave model.

Appendix B – Dimensionless definitions

References

All Tables

All Figures