Sage Journals: Discover world-class research

Abstract

Acid fracturing has been widely used as an industry practice in explored and developed carbonate reservoirs. It is very important to understand responses of reservoirs and improve production performance of a well due to the presence of fracture networks by stimulation treatments. Pressure transient analysis is one of the most effective diagnostic techniques available to enhance our understanding of natural and artificial-etched fracture behavior. This work presented a novel mathematical model for unsteady state flow of naturally fractured porous medium into multiple etched fractures intersecting a vertical well, and three different geometric shapes of matrix blocks containing slabs, cylinders and spheres were considered. The new solution was derived by using the Laplace transformation and the point source function integral approach. The polar coordinate transformation was used to deal with the radial distribution of arbitrary fracture number and angle. Then the model was validated by comparison with three published cases. Finally, type curves were plotted to identify flow regimes: linear flow, transitional flow, pseudoradial flow, and boundary dominant flow if the closed or constant pressure boundary exists. Furthermore, sensitivity analysis was investigated. The results showed that the acid-etched fracture parameters containing fracture number, fracture distribution and conductivity had a significant impact on pressure behavior at early times. However, natural fracture storativity coefficient and interporosity flow parameter mainly affected the transitional flow at intermediate times. Moreover, the shape of matrix blocks had a little influence on transient responses at intermediate times. It is found that multiple etched fractures existing near the wellbore consume less pressure drop and increase the productivity of a well as a whole.

Keywords

Carbonate reservoirs etched fractures vertical well pressure transient analysis unsteady state flow type curve

Introduction

Carbonate reservoirs possess about half of the world proven oil reserves and over 60% of world oil production, and have become increasingly important resources (Bagrintseva, 2015). The presence of natural fracture is typical feature of carbonate rocks, and thus affects many aspects of exploitation and production in fractured carbonate reservoirs, such as acid fracturing treatments, completion strategies, reservoir management and production (Wayne, 2008; Jiao et al., 2018). The spatial variability of fracture leads to a tremendous challenge to understand the behavior and performance of naturally fractured carbonate reservoirs when a vertical well with acid fracturing treatment is produced. Many scholars have proposed special approaches to better characterize flow behavior of naturally fractured carbonate formations. Analytical, semianalytical and numerical models have been presented to perform pressure transient analysis; at least two kinds of models have been developed to study the transient responses of such formations: the dual-porosity model and discrete fracture model.

The dual-porosity model that takes into account the fluid transfer between matrix and fractures includes two different conditions: the pseudosteady state and unsteady state flow. For these models, two key parameters, the storativity coefficient and interporosity flow parameter, are sufficient to characterize the behavior of the double medium. Since Warren and Root (1963) first proposed an idealized model that assumes a network of orthogonal, equally spaced feature of fracture system, with the pseudosteady-state flow between the matrix and fracture, more than a hundred authors have used the dual-porosity model under the pseudosteady state flow to investigate and interpret pressure transient responses of naturally fractured reservoirs in the petroleum engineering domain (Samaniego and Cinco-Ley, 2009). However, Kuchuk and Biryukov (2014) pointed out that the Warren and Root dual-porosity model is a fictitious homogeneous porous medium because it does not contain any fractures, and discovered that most well tests published in the literature do not exhibit the Warren and Root dual-porosity model behavior. On the other hand, the other dual-porosity model was widely used on the basis of the assumption of a transient exchange term between matrix and fracture. Several researchers have presented different geometric shapes of matrix blocks, such as slab, cylindrical and spherical matrix blocks, to describe the unsteady state exchange (de SWAAN, 1976; Kucuk and Sawyer 1980, Cinco-Ley and Samaniego 1982; Olarewaju and Lee 1986; Celis et al., 1994; Hassanzadeh et al., 2009; Babak and Azaiez 2014; Wang et al., 2015). The biggest difference between the pseudosteady state and nonsteady state flow is that the pressure derivative curve of the former has a deeper trough like a bell-shape than the latter one being relatively flat at the transition flow period.

In recent years, research interests have been growing to evaluate the naturally fractured reservoir on the use of the discrete fracture model. Several analytical or semianalytical solutions for understanding the flow behavior in discretely fractured porous media were presented by Izadi and Yildiz (2009), Zeng et al. (2012), Biryukov and Kuchuk (2012), Kuchuk and Biryukov (2014), Jia et al. (2015, 2016) and Chen et al. (2017). Their results showed that the pressure behavior of discretely fractured formations is considerably different from the classical model of Warren and Root’s pseudosteady state behavior. In some cases, the discrete fracture model exhibits a dual-porosity model behavior of unsteady state flow. Noushabadi et al. (2013) investigated the hydrodynamic response during a well test for uniform distribution of the fracture network though numerical simulation. All the above-mentioned studies of discretely fractured reservoirs require exact knowledge of fracture distribution, which results in difficultly establishing discrete fracture network models.

On the other hand, most producers from carbonate reservoirs are necessary for acidizing or acid fracturing to generate conductive flow paths connecting wellbore and reservoir because of extremely tight carbonate rocks. In an acidizing process, an acidic solution is injected to dissolve some of the carbonate rock, creating conductive channels called wormholes (Maheshwari and Balakotaiah, 2013). Liu et al. (2006) developed a wormhole model for cold heavy oil production. Wang et al. (2014) presented a vertical well intersected multiple wormholes in naturally fractured-vuggy porous media, and their model is based on multiple-continuum medium concept without considering the fracture conductivity. These wormhole models are not appropriate for multiple etched fractures. In an acid fracturing process, multiple etched non-smooth fracture surfaces leave open pathways, in addition to the wormholes created from the artificial fracture to the formation (Sarma et al., 2010). Deng et al. (2012) pointed out that the acid-etched fracture conductivity depends on permeability distribution, calcite percentage and rock property. But the mineralogy distribution dominates the etching patterns. Unlike hydraulic fractures with basically orthonormal and symmetrical distribution, high conductivity because of proppant placement, and longer fracture length, the etched fracture exhibits a radial distribution with shorter fracture length and lower conductivity due to fracture closure (Li et al., 2015; Lee et al., 2016; He et al., 2016; Qin et al., 2018). The presence of multiple etched fractures makes pressure transient analysis more complex and difficult, especially the radial distribution of multiple etched fractures. Generally, hydraulic fractures with orthogonality usually occur in a multifractured horizontal well, and the Green function and Newman product methods are often used to derive the solution. However, the polar coordinate transformation based on the point source function is very effective to solve the solution of multiple fractures intersecting a vertical well (Wang et al., 2014). A symmetrical single vertical fracture in double porosity reservoirs was studied by Cinco-Ley and Meng (1988). Liu et al. (2019) derived a mathematical model of a vertical fracture with varying conductivity. For multiple fractures connected to a well for acidizing or acid fracturing processes, the composite model with two zones is a popular way (Olarewaju and Lee, 1989; Bensalem and Tiab, 2007). Karasaki et al. (1988) to simplify the complicated problem of multiple fractures, and in the inner region the flow was assumed to be linear. However, in the outer region, the flow was assumed to be radial. In addition, the diagnostic plot of real-time monitoring of matrix acidizing was also proposed to calculate the evolving skin factor (Hill and Zhu, 1996). Recently, multiple fractures have been extensively explored using analytical approach by Fisher et al. (2005), Chen et al. (2014, 2017), Ren and Guo (2015) and Luo and Tang (2015, 2017). In numerical well test for multiple fractures, Choo and Wu (1987) studied the transient pressure behavior of multiple-fractured gas wells. Wan et al. (2016) used the finite element method to deal with fractures with branches.

Although considerable efforts have been dedicated to study either naturally fractured reservoirs or multiple fractures in homogeneous reservoirs, to the best of our knowledge, there is not a complete study that contains multiple fractures intersecting a vertical well in fractured carbonate reservoirs. The objective of our research is to seek a new solution for calculating pressure responses of a vertical well containing multiple fractures connected to the wellbore in carbonate reservoirs. Once the mathematical model is available based on the point source function, the next step is to solve the model, and plot type curves of pressures and pressure derivatives to identify the flow regimes of multiple fractures connected to a vertical well in carbonate reservoirs. Sensitive analysis is also discussed in detail on number of fractures, fracture distribution, fracture conductivity, storativity coefficient, interporosity flow parameter and different geometric shapes of matrix blocks.

Mathematical model and algorithm

Let us consider a vertical well intersected by multiple etched fractures with finite conductivity, which is produced in a carbonate reservoir as indicated in Figure 1. The dual-porosity model under the unsteady state flow is adopted to describe the fluid exchange between the matrix and fracture system for naturally fractured formations. Three different geometric shapes of matrix blocks are assumed: linear flow for slab shape, radial flow for cylindrical and spherical shape as shown in Figure 1(a). The other assumptions are made as follows:

Figure 1.

Schematic diagram. (a) Actual fractured reservoirs idealized as slab-shaped, cylindrical-shaped and spherical-shaped. (b) The system of etched fractures, natural fractures and matrix. (c) Flow model of the system.

All rock properties including permeability and porosity are constant. Fluid in the formation is slightly compressible. The formation thickness is constant.

Only oil flow is assumed as single phase at isothermal condition and Darcy flow.

Three types of matrix blocks are considered: slab-shape, cylindrical-shape and spherical-shape. The matrix blocks act as a uniformly distributed source in the fractured medium.

The number of etched fractures is arbitrary. Fluids first exchange from the matrix system into the natural fracture system, then enter into multiple etched fractues, and finally go into the well, as shown in Figure 1(c). The total rate of multiple etched fractures is equal to the total rate of the production well.

The external boundaries of the reservoirs are infinite, closed or constant pressure.

Point source function

According to the Green’s and source functions theory, the point source solution is the basic of the line source solution (Gringarten et al., 1974). We first should establish a point source function of naturally fractured carbonate reservoirs. The use of these dimensionless variables allows us to present the solution for the flow equations in a simple general form. Appendix 1 gives these dimensionless definitions in detail.

Flow model in reservoir

The partial differential equation of pressure diffusion for matrix with the assumption of an unsteady state flow in dimensionless form can be written as follows (Hassanzadeh et al., 2009; Babak and Azaiez 2014; Wang et al., 2015)

For slab-shaped matrix blocks

\frac{\partial^{2} p_{m D}}{\partial x_{D}^{2}} = \frac{3 (1 - ω)}{λ} \frac{\partial p_{m D}}{\partial t_{D}}

(1)

For cylindrical-shaped matrix blocks

\frac{\partial^{2} p_{m D}}{\partial r_{D}^{2}} + \frac{1}{r_{D}} \frac{\partial p_{m D}}{\partial r_{D}^{}} = \frac{8 (1 - ω)}{λ} \frac{\partial p_{m D}}{\partial t_{D}^{}}

(2)

For spherical-shaped matrix blocks

\frac{\partial^{2} p_{m D}}{\partial r_{D}^{2}} + \frac{2}{r_{D}} \frac{\partial p_{m D}}{\partial r_{D}^{}} = \frac{15 (1 - ω)}{λ} \frac{\partial p_{m D}}{\partial t_{D}^{}}

(3)

The initial and boundary conditions for the matrix media are given by

p_{m D} (x_{D}, t_{D} = 0) = 0

(4)

\frac{\partial p_{m D} (x_{D} = 0, t_{D})}{\partial x_{D}} = 0

(5)

p_{m D} (x_{D} = 1, t_{D}) = p_{f D}

(6)

where

x_{D} = x / (h_{m} / 2)

and

λ = 12 λ^{*}

(slab shape);

r_{D} = r / R_{m}

and

λ = 8 λ^{*}

(cylindrical shape);

r_{D} = r / R_{m}

and

λ = 15 λ^{*}

(spherical shape). For radial flow of cylindrical shape and spherical shape, using r_D substituting for x_D.

Applying the Laplace transformation to equations (1) to (3), we obtain the following system of equations

\frac{d^{2} {\tilde{p}}_{m D}}{d x_{D}^{2}} = \frac{3 (1 - ω)}{λ} s {\tilde{p}}_{m D}

(7)

\frac{d^{2} {\tilde{p}}_{m D}}{d r_{D}^{2}} + \frac{1}{r_{D}} \frac{d {\tilde{p}}_{m D}}{d r_{D}^{}} = \frac{8 (1 - ω)}{λ} s {\tilde{p}}_{m D}

(8)

\frac{d^{2} {\tilde{p}}_{m D}}{d r_{D}^{2}} + \frac{1}{r_{D}} \frac{d {\tilde{p}}_{m D}}{d r_{D}^{}} = \frac{15 (1 - ω)}{λ} s {\tilde{p}}_{m D}

(9)

The above three ordinary-differential equations can be solved using the inner and outer boundary conditions. However, we only need to acquire pressure gradients within the matrix blocks, and evaluate the flux term at the matrix-fracture interface, which is accounting for the unsteady state flow. These solutions of pressure gradients in the Laplace space are given by

{(\frac{d {\tilde{p}}_{m D}}{d x_{D}^{}})}_{x_{D} = 1} = \sqrt{\frac{3 (1 - ω) s}{λ}} \tanh \sqrt{\frac{3 (1 - ω) s}{λ}} \cdot {\tilde{p}}_{f D}

(10)

{(\frac{d {\tilde{p}}_{m D}}{{dr}_{D}^{}})}_{r_{D} = 1} = \frac{\sqrt{\frac{8 (1 - ω) s}{λ}} I_{1} (\sqrt{\frac{8 (1 - ω) s}{λ}})}{I_{0} (\sqrt{\frac{8 (1 - ω) s}{λ}})} {\tilde{p}}_{f D}

(11)

{(\frac{d {\tilde{p}}_{m D}}{{dr}_{D}^{}})}_{r_{D} = 1} = (\sqrt{\frac{15 (1 - ω) s}{λ}} \coth \sqrt{\frac{15 (1 - ω) s}{λ}} - 1) {\tilde{p}}_{f D}

(12)

It is noted that the pressure gradient term of matrix blocks contains the fracture pressure that can be used to couple the fracture equation.

Flow model in fracture

The governing equation for fracture in dimensionless form can be written as

For slab-shaped matrix blocks

\frac{\partial^{2} p_{f D}}{\partial r_{D}^{2}} + \frac{1}{r_{D}} \frac{\partial p_{f D}}{\partial r_{D}^{}} - \frac{λ}{3} {(\frac{\partial p_{m D}}{\partial x_{D}})}_{x_{D} = 1} = ω \frac{\partial p_{f D}}{\partial t_{D}^{}}

(13)

For cylindrical-shaped matrix blocks

\frac{\partial^{2} p_{f D}}{\partial r_{D}^{2}} + \frac{1}{r_{D}} \frac{\partial p_{f D}}{\partial r_{D}^{}} - \frac{λ}{4} {(\frac{\partial p_{m D}}{\partial r_{D}})}_{r_{D} = 1} = ω \frac{\partial p_{f D}}{\partial t_{D}^{}}

(14)

For spherical-shaped matrix blocks

\frac{\partial^{2} p_{f D}}{\partial r_{D}^{2}} + \frac{1}{r_{D}} \frac{\partial p_{f D}}{\partial r_{D}^{}} - \frac{λ}{5} {(\frac{\partial p_{m D}}{\partial r_{D}})}_{r_{D} = 1} = ω \frac{\partial p_{f D}}{\partial t_{D}^{}}

(15)

Equations (13) to (15) can be converted into the Laplace domain as given by

\frac{d^{2} {\tilde{p}}_{f D}}{d r_{D}^{2}} + \frac{1}{r_{D}} \frac{d {\tilde{p}}_{f D}}{d r_{D}^{}} - \frac{λ}{3} {(\frac{d {\tilde{p}}_{m D}}{d x_{D}})}_{x_{D} = 1} = ω s {\tilde{p}}_{f D}

(16)

\frac{d^{2} {\tilde{p}}_{f D}}{d r_{D}^{2}} + \frac{1}{r_{D}} \frac{d {\tilde{p}}_{f D}}{d r_{D}^{}} - \frac{λ}{4} {(\frac{d {\tilde{p}}_{m D}}{d r_{D}})}_{r_{D} = 1} = ω s {\tilde{p}}_{f D}

(17)

\frac{d^{2} {\tilde{p}}_{f D}}{d r_{D}^{2}} + \frac{1}{r_{D}} \frac{d {\tilde{p}}_{f D}}{d r_{D}^{}} - \frac{λ}{5} {(\frac{d {\tilde{p}}_{m D}}{d r_{D}})}_{r_{D} = 1} = ω s {\tilde{p}}_{f D}

(18)

Substituting equation (10) into equation (16), equation (11) into equation (17) and equation (12) into equation (18), respectively, three equations of different shapes become a general form

\frac{d^{2} {\tilde{p}}_{f D}}{d r_{D}^{2}} + \frac{1}{r_{D}} \frac{d {\tilde{p}}_{f D}}{d r_{D}^{}} = s f (s) {\tilde{p}}_{f D}

(19)

where

f (s) = {\begin{matrix} \frac{λ}{3 s} \sqrt{\frac{3 (1 - ω) s}{λ}} \tanh \sqrt{\frac{3 (1 - ω) s}{λ}} + ω, slab shape \\ \frac{λ}{4 s} \sqrt{\frac{8 (1 - ω) s}{λ}} I_{1} (\sqrt{\frac{8 (1 - ω) s}{λ}}) / I_{0} (\sqrt{\frac{8 (1 - ω) s}{λ}}) + ω, cylindrical shape \\ \frac{λ}{5 s} (\sqrt{\frac{15 (1 - ω) s}{λ}} \coth \sqrt{\frac{15 (1 - ω) s}{λ}} - 1) + ω, spherical shape \end{matrix}

(20)

where I₀ and I₁ are modified Bessel functions of the first kind of order zero and one, respectively.

The initial and boundary conditions in the Laplace space for fractured media are given by

{\tilde{p}}_{f D} (x_{D}, s = 0) = 0

(21)

\frac{d {\tilde{p}}_{f D} (x_{D} = 0, s)}{d x_{D}} = - \frac{1}{s}

(22)

{\tilde{p}}_{f D} (x_{D} \to \infty, s) = 0

(23)

Let the polar coordinates of the any point (x, y) and the source point (xʹ, yʹ) be given, respectively, by ( $r$ , $θ$ ) and ( ${r^{'}}_{}$ , $θ'$ ). The general solution of equation (19) can be solved with the above initial and boundary conditions given

{\tilde{p}}_{f D} = {\tilde{q}}_{D} (s) K_{0} (ξ_{D} \sqrt{s f (s)})

(24)

where

ξ_{D} = \sqrt{r_{D}^{2} + {r^{'}}_{D}^{2} - 2 r_{D} {r^{'}}_{D} \cos (θ - θ')}

(25)

where K₀ is the modified Bessel functions of the second kind of order zero.

Equation (24) represents a point source function in the Laplace domain.

Single fracture solution

Let us redefine the dimensionless variables $r_{D} = r / x_{f}$ and ${r^{'}}_{D} = r^{'} / x_{f}$ . According to the point source function and the principle of superposition (Ozkan and Raghavan, 1991), integrating the point source function equation (24) along the fracture plane leads to a line solution of the single vertical fracture

{\tilde{p}}_{f D} = \int_{0}^{1} {\tilde{q}}_{D} ({r^{'}}_{D}, s) K_{0} (ξ_{D} \sqrt{s f (s)}) d {r^{'}}_{D}

(26)

Assuming that the flux is uniformly distributed along the fracture, then

{\tilde{p}}_{f D} = {\tilde{q}}_{D} (s) \int_{0}^{1} K_{0} (ξ_{D} \sqrt{s f (s)}) d {r^{'}}_{D}

(27)

In order to obtain the solution of single infinite conductivity fracture, ignoring the pressure drop along the fracture, the superposition integral can be utilized to approximate the infinite conductivity solution by integrating along the fracture plane using the uniform flux solution (Craig and Blasingme, 2006; Restrepo and Tiab, 2009)

{\tilde{p}}_{f D}^{\inf} = \int_{0}^{1} {\tilde{q}}_{D} (s) \int_{0}^{1} K_{0} (ξ_{D} \sqrt{s f (s)}) d r_{D}^{'} d r_{D}

(28)

In fact, the etched fracture conductivity is generally finite and the pressure drop along the fracture is varying. For a single finite conductivity fracture, Wilkinson (1989) and Riley et al. (1991) proposed the conductivity influence function concept that could be combined into the infinite conductivity solution. The conductivity influence function is similar to a dynamic skin representing the impact of finite conductivity. Later, Wang et al. (2014), Wang et al. (2017) and Xu et al. (2017) applied the influence function to other aspects of multiple fracture studies. Recently, Luo et al. (2018) defined directly the difference of dimensionless wellbore pressure between finite conductivity fracture and infinite conductivity fracture as conductivity influence function, and derived the mathematical expression in detail. Thus, we adopt the rapid and accurate method of influence function to consider finite conductivity fracture, avoiding some complicated computations like Cinco-Ley and Meng (1988) and Riley et al. (1991). The dimensionless pressure distribution in the Laplace domain can be expressed as

s {\tilde{p}}_{f D} (r_{D}, s) = s {\tilde{p}}_{_{f D}}^{inf} (r_{D}, s) + s \tilde{f} (C_{f D})

(29)

where the influence function

\tilde{f} (C_{f D})

in this work may be written as

s \tilde{f} (C_{f D}) = 2 π \sum_{n = 1}^{\infty} \frac{1}{n^{2} π^{2} C_{f D} + 2 \sqrt{n^{2} π^{2} + s}} + \frac{0.4063 π}{π (C_{f D} + 0.8997) + 1.6252 s}

(30)

Therefore, a single finite conductivity fracture solution can be obtained by combining infinite conductivity solution equation (28) with the influence function equation (30) in equation (29) as

s {\tilde{p}}_{f D} = s \int_{0}^{1} {\tilde{q}}_{D} (s) \int_{0}^{1} K_{0} (\sqrt{r_{D}^{2} + {r^{'}}_{D}^{2} - 2 r_{D} {r^{'}}_{D} \cos (θ - θ')} \sqrt{s f (s)}) d r_{D}^{'} d r_{D} + s \tilde{f} (C_{f D})

(31)

Multiple fractures solution

Considering n multiple etched fractures as shown in Figure 1(b), the principle of superposition (Wang et al., 2014) is used to calculate the pressure response caused by the etched fractures due to the linear partial differential equation on pressure. Thus, on the basis of equation (31), we can compute the dimensionless pressure at any point of the reservoir in the Laplace domain due to the n etched fractures producing

s {\tilde{p}}_{f D} = \sum_{j = 1}^{n} s [\int_{0}^{1} {\tilde{q}}_{D j} (s) \int_{0}^{1} K_{0} (ξ_{D j} \sqrt{s f (s)}) d r_{D}^{'} d r_{D}] + δ_{j} \tilde{f} (C_{f D})

(32)

where

ξ_{D j} = \sqrt{r_{D}^{2} + {r^{'}}_{D}^{2} - 2 r_{D} r_{D}^{'} \cos (θ - θ_{j})}

(33)

Using equation (32), the dimensionless pressure drop of the kth etched fracture can be obtained by

s {\tilde{p}}_{fDk} = \sum_{j = 1}^{n} s [\int_{0}^{1} {\tilde{q}}_{D j} (s) \int_{0}^{1} K_{0} (ξ_{Djk} \sqrt{s f (s)}) d r_{D}^{'} d r_{D}] + δ_{j} \tilde{f} (C_{f D}) k = 1, 2, \dots, n

(34)

where

ξ_{Djk} = \sqrt{r_{D}^{2} + {r^{'}}_{D}^{2} - 2 r_{D} r_{D}^{'} \cos (θ_{k} - θ_{j})}

(35)

and

δ_{j} = {\begin{matrix} 1, & j \neq k \\ 0, & j = k \end{matrix}

(36)

Because the sum of all the etched fracture rates should be equal to the total flow rate of the vertical well, we have

\sum_{j = 1}^{n} s {\tilde{q}}_{fDj} (s) = 1

(37)

Applying equation (34) to all etched fractures, the assumption of the wellbore pressure requires equivalent pressure at the location where fluid enters the wellbore from each fracture. Thus, a n + 1-order matrix is obtained as follows

[\begin{matrix} {\tilde{p}}_{w D 11} & {\tilde{p}}_{w D 12} & \dots & {\tilde{p}}_{w D 1 n} & - 1 \\ {\tilde{p}}_{w D 21} & {\tilde{p}}_{w D 22} & \dots & {\tilde{p}}_{w D 2 n} & - 1 \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ {\tilde{p}}_{wDn 1} & {\tilde{p}}_{wDn 2} & \dots & {\tilde{p}}_{wDnn} & - 1 \\ 1 & 1 & \dots & 1 & 0 \end{matrix}] [\begin{matrix} {\tilde{q}}_{f D 1} \\ {\tilde{q}}_{f D 2} \\ ⋮ \\ {\tilde{q}}_{fDn} \\ {\tilde{p}}_{w D} \end{matrix}] = [\begin{matrix} 0 \\ 0 \\ ⋮ \\ 0 \\ 1 / s \end{matrix}]

(38)

The n + 1 unknown variables can be solved by applying Gauss elimination. Then the dimensionless wellbore pressure ${\tilde{p}}_{w D}$ in the Laplace domain may be obtained. Finally, we can get the dimensionless transient pressure $p_{w D}$ in the time domain by using the Stehfest numerical algorithm (Stehfest, 1970). At the same time, the dimensionless pressure derivative was computed for creating log–log diagnostic plots and examining transient flow regimes.

Model validation

In this section, we will compare the results of our model with the results from three literature data published. Two cases of the fracture were considered: infinite conductivity and finite conductivity. One example is that Gringarten et al. (1974) derived a single infinite conductivity vertical fracture solution and provided detailed data of p_wD versus t_D. The second example is that Ozkan and Raghavan (1991) provided a library of uniform flux solutions in the Laplace domain for different wellbore configurations and boundary conditions. The infinite conductivity solution can be obtained by evaluating the uniform flux solution at an equivalent average pressure point (x_D=0.732). Figure 2(a) shows the comparison results of infinite conductivity fracture in this paper with the results from Gringarten et al. and Ozkan and Raghavan. We note excellent agreement for the whole trend. The advantage of our new model is that it is able to deal with more than two fractures. The third example is that Cinco-Ley et al. (1976) developed a finite conductivity fracture solution and provided output data under different conductivity values. Figure 2(b) illustrates the wellbore pressure comparison of this work and the Cinco-Ley method for finite conductivity case, and we also observe excellent agreement for all curve trends. Whatever infinite or finite conductivity, our results are consistent with the data in the literature. All of the three verifications also demonstrate that the finite conductivity solution by adding the influence function into the infinite conductivity solution is accurate and can be used to analyze flow behavior.

Figure 2.

(a) Verification of wellbore pressure under infinite conductivity for data of Gringarten et al. (1974) and Ozkan and Raghavan (1991) and (b) Verification of wellbore pressure under different finite conductivity conditions for data of Cinco (1978).

Results and discussion

In this section, flow regimes of type curves first are analyzed. In order to distinguish the effect of the outer boundary, three cases are considered: infinite, closed and constant pressure. Then, we will discuss in detail the effect of sensitive parameters on the pressure and pressure derivative response. These key parameters include the number of fractures, fracture distribution, fracture conductivity, fracture storativity coefficient, interporosity flow parameter and different geometric configuration of matrix blocks. The changing characteristics of pressure transient behavior for a vertical well with multiple etched fractures in carbonate reservoirs are investigated as follows.

Flow regimes

The input parameters are necessary to generate these types of curves. As a base case, these values are each fracture conductivity C_fD=10π, fracture storativity coefficient ω = 0.01, interporosity flow parameter λ = 1 × 10⁻⁴ and dimensionless reservoir radius r_eD=3000, with the assumption of three equal fractures, an angle θ of 120° between arbitrary two fractures and slab-shape matrix blocks. The type curves of the vertical well intersecting multiple fractures in this paper are plotted as shown in Figure 3. Obviously, four flow regimes would be identified from dimensionless pressure derivative curves: linear flow at an early time, transitional flow at an intermediate time, pseudoradial flow at a late time, and probable boundary dominant flow if the presence of closed or constant pressure boundaries exists. The first flow regime is multiple fracture effect near the vertical well characterized by a straight line with unit slope. Note that this is different from wellbore storage effect. The second regime is the transitional flow period with a relatively flat trough, of which the duration depends mainly on the double media nature of carbonate reservoirs. The third flow regime is the pseudoradial flow period as the fluid converges radially to the entire system, which exhibits a horizontal line of 0.5 on the pressure derivative curve. The fourth flow regime is the boundary dominant flow when the pressure derivative goes upward for closed outer boundary and goes downward for constant pressure boundary.

Figure 3.

Type curves with different outer boundaries.

Effect of number of fracture

It is assumed that the conductivity of each fracture is infinite, the shape of matrix blocks is slab, multiple fractures are symmetrical distributions and the outer boundary is infinite. Other parameters (ω and λ) are constant. We design seven cases of fracture numbers equal to 2, 3, 4, 5, 6, 7 and 8, and their corresponding angles are 90°, 120°, 90°, 72°, 60°, 51.4°and 45°, respectively. The effect of fracture number on pressure and pressure derivative curves is plotted in Figure 4. As can be seen, etched fractures only influence at early time period and the pressure behavoir is a function of the fracture number. The impact would last until the start of the intermediate time transitional flow. The dimensionless pressure and pressure derivative at early times decrease as the fracture number increases, which indicates that the more fracture number, the same production rate needs consume less energy. Furthermore, the difference in the wellbore pressure response decreases as the fracture number exceeds four. That means that the interference between etched fractues occurs. It is noted that the more fracture number generates the higher conductivity of the articalfical fracture system near the wellbore as a whole.

Figure 4.

(a) Effect of fracture number on the pressure and (b) Effect of fracture number on the pressure derivative.

As expected, etched fractues formed during the process of carbonate acidizing can transform the formation radial flow into the linear flow. That is to say, artifical fractures connected the wellbore improve flow capability near the well, and the well productivity would be enhanced by multiple etched fractues under the same production pressure drop.

Effect of fracture distribution

In Figure 5, two fractures of various angles are considered. One fracture along the direction of the x axis is assumed and the other fracture keeps an arbitrary angle θ with the x direction. Five different cases of fracture angles equal to 1°, 5°, 15°, 60° and 180°, respectively, are taken into consideration to analyze their effects on pressure behavior. As shown in Figure 5(b), the expected three flow regimes containing the linear flow, the transient state crossflow and the pseudoradial flow period are observed in all the cases. Similar to the effect of fracture numbers, the pressure and pressure derivative curves are affected at early times by the angle. As the angle increases, the pressure drop decreases, and the duration of linear flow ends later. Moreover, at the first half part of early times (1 × 10⁻⁶<t_D<1 × 10⁻⁴), the incremental value of pressure transient response reduces firstly rapidly then slowly, and at the second half part of early times (1 × 10^-<t_D<1 × 10⁻²), the incremental value of pressure and pressure derivative decreases firstly slowly then swiftly. That is the reason that the interference between two fractures becomes weaker with the increase of the angle.

Figure 5.

(a) Effect of the angle of two etched fractures on the pressure and (b) Effect of the angle of two etched fractures on the pressure derivative.

To further illustrate the effect of fracture distribution on pressure behavior, three different scenarios of etched fracture distributions are considered when the fracture number keeps the same. With the assumption of four symmetrical fractures, we design three patterns of the fracture distribution with 30°, 45°, 90° of the minimum angle. Figure 6 demonstrates that the patterns of fracture distribution only influence the early regime. As the minimum angle increases, the pressure derivative moves slightly down. Because the interference between these fractures becomes weak with the increase of the minimum angle. On the other hand, why the influence is so slight. We have recognized that the difference in the pressure derivative response is relatively small when the fracture number is beyond four according to the above-mentioned effect of fracture number. That is to say, four fractures are enough for the acid treatment around the vertical well when the minimum angle exceeds certain value.

Figure 6.

Effect of fracture distribution on the pressure derivative for four etched fractures.

Effect of fracture conductivity

Figure 7 shows the effect of fracture conductivity on the pressure and pressure derivative. For three etched fractures symmetrically with the angle of 120°, linear flow can be identified when C_fD>π. Similarly, fracture conductivity only affects the early period. As fracture conductivity increases, the pressure drop and pressure derivative curves decrease. The pressure derivative response for C_fD=π and C_fD=10π is close to each other at early times. However, the derivative behavior for C_fD=0.1π is obviously different from the other. The pressure derivative has an upward trend at late time because of closed outer boundary when we change the infinite condition into the closed condition.

Figure 7.

Effect of fracture conductivity on the pressure and pressure derivative for three etched fractures.

Effect of fracture storativity coefficient

The storativity coefficient of the fracture system represents the relative capacity of fluid stored in a fracture system. According to the definitions of storativity coefficients of fractures and matrix, their sum is equal to one. Figure 8 shows the effect of fracture storativity coefficient on pressure and pressure derivative responses in naturally fractured carbonate reservoirs. The maximum fracture storativity coefficient is 100 times higher than the minimum one. It is found that the ω not only decides the duration and depth of the transitional flow period between matrix blocks and natural fractures, but also has a significant effect on the starting time of the early flow period that is different from the early wellbore storage flow. It can be seen that a smaller ω will lead to wider and deeper transition flow regimes. The concave segment of the transition flow reflects the characteristics of naturally fractured carbonate reservoirs under the nonsteady state. For ω = 0, the type curve becomes the same as a homogeneous reservoir.

Figure 8.

Effect of fracture storativity coefficient on the pressure and pressure derivative for three etched fractures (C_fD=10π).

Effect of interporosity flow parameter

The interporosity flow parameter between a matrix system and fracture system represents the fluid starting time of the matrix system to the fracture system. The λ values are set as 1 × 10⁻³, 1 × 10⁻⁴ and 1 × 10⁻⁵, respectively, while keeping the other parameters constant. The wellbore pressure and pressure derivative responses for three cases are displayed in Figure 9. As the λ increases, the pressure derivative uplifts earlier at intermediate times, then keeps constant for a period of time, which implies that the λ values mainly influence the starting time of nonsteady state crossflowing from matrix to fractures.

Figure 9.

Effect of interporosity flow parameter on the pressure and pressure derivative for three etched fractures (C_fD=10π).

Effect of different geometric shape of matrix blocks

Figure 10 shows a log–log graph of dimensionless pressure and pressure derivative versus dimensionless time for three different geometric shapes of matrix blocks. At early times, the pressure responses are independent of geometric shapes of matrix blocks, and exhibit linear flow, recognizable by unit slope. At intermediate times, all curves almost overlay except being a little diverging. The slab shape has the minimum point on the trough while the cylinder and sphere shape are almost similar. The reason can be explained that the slab shape is assumed to be linear flow of low resistance, while the cylindrical and spherical shapes are assumed to be radial flow.

Figure 10.

Effect of geometric configuration of matrix blocks on the pressure and pressure derivative for three etched fractures (C_fD=10π).

Conclusions

The mathematical model of multiple etched fractures connected a vertical well is developed in naturally fractured carbonate reservoirs. An unsteady-state exchange model is applied to describe fluid exchange between the matrix and fracture system. Three types of matrix blocks containing slab, cylinder and sphere are assumed. The proposed model is an extension of the classic dual-porosity model based on a continuum medium concept.

New analytical solutions are obtained by using Laplace transformation, point source integral and Stehfest inversion methods. The fracture conductivity function is added to the expression of infinite conductivity solution and thus can easily deal with various conductivity conditions.

Type curves were generated to identify the flow characteristics. Four flow regimes may occur: linear flow at early times, transitional flow at intermediate times, pseudoradial flow at late times, and perhaps boundary dominant flow at late times if the closed or constant pressure boundaries exist.

The sensitive parameters, such as fracture number, fracture distribution, fracture conductivity, fracture storativity coefficient, interporosity flow parameter and different geometric shapes of matrix blocks, are investigated. The former three represent the effect of artificial fractures on transient pressure near the wellbore. Thus, pressure and pressure derivative curves only change at early times. The natural fracture parameters of the storativity coefficient and interporosity flow parameter are responsible for the intermediate time period in naturally fractured medium. The geometric shapes of matrix blocks including slabs, cylinders and spheres have a weak effect on pressure behavior.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors greatly appreciate the financial support from the Important National Science and Technology Specific Projects of China (Grant No. 2017ZX05030-002).

References

Babak

Azaiez

(2014) Unified fractional differential approach for transient interporosity flow in naturally fractured media. Advances in Water Resources 74: 302–317.

Bagrintseva

(2015) Carbonate Reservoir Rocks. Hoboken, NJ: Wiley.

Bensalem

Tiab

(2007) Interpretation of post deep-acidizing pressure tests in a naturally fractured reservoir. Paper 107228 presented at the 2007 SPE Latin American and Caribbean Petroleum Engineering Conference, Buenos Aires, Argentina, 15–18 April.

Biryukov

Kuchuk

(2012) Transient pressure behavior of reservoirs with discrete conductive faults and fractures. Transport in Porous Media 95: 239–268.

Celis

Silva

Ramones

, et al. (1994) A new model for pressure transient analysis in stress sensitive naturally fractured reservoirs. SPE Advanced Technology Series 2(1): 126–135.

Chen

Liao

Huang

, et al. (2014) Productivity estimations for vertically fractured wells with asymmetrical multiple fractures. Journal of Natural Gas Science and Engineering 21: 1048–1060.

Chen

Liao

Sepehrnoori

, et al. (2017) A semi-analytical model for pressure transient analysis of fractured wells in unconventional plays with arbitrarily distributed fracture networks. Paper 187290 presented at the SPE annual technical conference and exhibition, San Antonio, Texas, 9–11 October.

Chen

Liao

Zhao

, et al. (2017) A comprehensive productivity equation for multiple fractured vertical wells with non-linear effects under steady-state flow. Journal of Petroleum Science and Engineering 149: 9–24.

Choo

(1987) Transient pressure behavior of multiple-fractured gas wells. Paper 16398 presented at the SPE/DOE low permeability reservoir symposium, Denver, Colorado, 18–19 May.

10.

Cinco-Ley

Meng

(1988) Pressure transient analysis of wells finite conductivity vertical fractures in double porosity reservoirs. Paper 18172 presented at the 63rd annual technical conference and exhibition, Houston, TX, 2–5 October.

11.

Cinco-Ley

Samaniego

(1982) Pressure transient analysis for naturally fractured reservoirs. Paper 11026 presented at the 57th annual fall technical conference and exhibition, New Orleans, LA, 26–29 September.

12.

Cinco-Ley

Samaniego

Dominguez

(1976) Transient pressure behavior for a well with a finite-conductivity vertical fracture. Society of Petroleum Engineers Journal 18(4): 253–264.

13.

Craig

Blasingme

(2006) Constant-rate drawdown solutions derived for multiple arbitrarily oriented uniform-flux, infinite-conductivity, or finite-conductivity fractures in an infinite-slab reservoir. Paper 100586 presented at the 2006 SPE gas technology symposium, Calgary, Alberta, Canada, 15–17 May.

14.

Deng

Mou

Hill

, et al. (2012) A new correlation of acid-fracture conductivity subject to closure stress. SPE Production & Operations 27(2): 158–169.

15.

de Swaan OA (1976) Analytic solutions for determining naturally fractured reservoir properties by well testing. Society of Petroleum Engineers Journal 16(3): 117–122.

16.

Fisher

Wright

Davidson

, et al. (2005) Integrating fracture-mapping technologies to improve stimulations in the Barnett shale. SPE Production & Facilities 20(2): 85–93.

17.

Gringarten

Ramey

Raghavan

(1974) Unsteady-state pressure distributions created by a well with a single infinite-conductivity vertical fracture. Society of Petroleum Engineers Journal 14(4): 347–360.

18.

Hassanzadeh

Pooladi-Darvish

Atabay

(2009) Shape factor in the drawdown solution for well testing of dual-porosity systems. Advances in Water Resources 32(11): 1652–1663.

19.

Cheng

, et al. (2016) A semianalytical methodology to diagnose the locations of underperforming hydraulic fractures through pressure-transient analysis in tight gas reservoir. SPE Journal 22(3): 1–16.

20.

Hill

Zhu

(1996) Real-time monitoring of matrix acidizing including the effects of diverting agents. SPE Production & Facilities 11(2): 95–101.

21.

Izadi

Yildiz

(2009) Transient flow in discretely fractured porous media. SPE Journal 14(2): 362–373.

22.

Jia

Cheng

Huang

, et al. (2015) Transient behavior of complex fracture networks. Journal of Petroleum Science and Engineering 132: 1–17.

23.

Jia

Cheng

Huang

, et al. (2016) A semi-analytical model for the flow behavior of naturally fractured formations with multi-scale fracture networks. Journal of Hydrology 537: 208–220.

24.

Jiao

, et al. (2018) Study on the effects of fracture on permeability with pore-fracture network model. Energy Exploration & Exploitation 36(6): 1556–1565.

25.

Karasaki

Long

Witherspoon

(1988) A new analytic model for fracture-dominated reservoirs. SPE Formation Evaluation 3(1): 242–250.

26.

Kuchuk

Biryukov

(2014) Pressure-transient behavior of continuously and discretely fractured reservoirs. SPE Reservoir Evaluation & Engineering 17(1): 82–97.

27.

Kucuk

Sawyer

(1980) Transient flow in naturally fractured reservoirs and its application to Devonian gas shales. Paper 9397 presented at the 55th annual fall technical conference and exhibition, Dallas, Texas, 21–24 September.

28.

Lee

Park

Shin

, et al. (2016) Efficient production estimation for a hydraulic fractured well considering fracture closure and proppant placement effects. Energy Exploration & Exploitation 34(4): 643–658.

29.

Jiao

Yue

, et al. (2015) Natural fractures and their effects on reservoir reconstruction in Lower Cambrian Shale, southeast Chongqing, China. Energy Exploration & Exploitation 33(6): 769–784.

30.

Liu

, et al. (2019) A mathematical model and semi-analytical solution for transient pressure of vertical fracture with varying conductivity in three crossflow rectangular layer. Energy Exploration & Exploitation 37(1): 230–250.

31.

Liu

Zhao

Jin

(2006) Coupled reservoir/wormholes model for cold heavy oil production wells. Journal of Petroleum Science and Engineering 50: 258–268.

32.

Luo

Tang

(2015) Pressure-transient analysis of multiwing fractures connected to a vertical wellbore. SPE Journal 20(2): 360–367.

33.

Luo

Wang

Tang

, et al. (2017) Productivity of multiple fractures in a closed rectangular reservoir. Journal of Petroleum Science and Engineering 157: 232–247.

34.

Luo

Wang

Liu

, et al. (2018) A simple and accurate calculation method for finite-conductivity fracture. Journal of Petroleum Science and Engineering 161: 590–598.

35.

Maheshwari

Balakotaiah

(2013) Comparison of carbonate HCL acidizing experiments with 3D simulations. SPE Production & Operations 28(4): 402–413.

36.

Noushabadi

Jourde

Massonnat

(2013) Evaluation of the inter-porosity exchange coefficient inferred from well test in a fractured reservoir with a non-uniform distribution pattern: A numerical study. Journal of Petroleum Science and Engineering 108: 82–95.

37.

Olarewaju

Lee

(1986) New pressure transient analysis model for dual-porosity reservoirs. Paper 15634 presented at the 61st annual technical conference and exhibition, New Orleans, LA, 5–8 October.

38.

Olarewaju

Lee

(1989) A comprehensive application of a composite reservoir model to pressure-transient analysis. SPE Reservoir Engineering 4(3): 325–331.

39.

Ozkan

Raghavan

(1991) New solutions for well-test-analysis problem: Part 1-analytical considerations. SPE Formation Evaluation 6(3): 359–3668.

40.

Qin

Cheng

, et al. (2018) An innovative model to evaluate fracture closure of multi-fractured horizontal well in tight gas reservoir based on bottom-hole pressure. Journal of Natural Gas Science and Engineering 57: 295–304.

41.

Ren

Guo

(2015) Performance of vertical fractured wells with multiple finite-conductivity fractures. Journal of Geophysics and Engineering 12: 978–987.

42.

Restrepo

Tiab

(2009) Multiple fractures transient response. Paper 121594 presented at the 2009 SPE Latin American and Caribbean petroleum engineering conference, Cartagena, Colombia, 31 May–3 June.

43.

Riley

Brigham

Horne

(1991) Analytic solutions for elliptical finite-conductivity fractures. Paper 22656 presented at the 66th annual technical conference and exhibition, Dallas, TX, USA, 6–9 October.

44.

Samaniego

Cinco-Ley

(2009) Transient well testing. In: Kamal M (ed) Naturally Fractured Reservoirs. Chap. 10. Richardson, Texas: Monograph Series, SPE, pp. 221–280.

45.

Sarma

Roy

Roul

, et al. (2010) Effect of acid type on initial response and sustenance of production in acid fracturing treatment: A field treatment evaluation. Paper 127339 presented at the SPE oil and gas India conference and exhibition, Mumbai, India, 20–22 January.

46.

Stehfest

(1970) Numerical inversion of Laplace transform. Communications of the ACM 13: 47–49.

47.

Wan

Liu

, et al. (2016) A numerical approach for pressure transient analysis of a vertical well with complex fractures. Acta Mechanica Sinica 32(4): 640–648.

48.

Wang

Dong

(2017) Rate decline curves analysis of multiple-fractured horizontal wells in heterogeneous reservoirs. Journal of Hydrology 553: 527–539.

49.

Wang

Luo

, et al. (2014) Analytical modeling of flow behavior for wormholes in naturally fractured-vuggy porous media. Transport in Porous Media 105(3): 539–558.

50.

Wang

Ding

, et al. (2015) Approximate analytical-pressure studies on dual-porosity reservoirs with stress-sensitive permeability. SPE Reservoir Evaluation & Engineering 18(4): 523–533.

51.

Wang

Luo

Hou

, et al. (2014) Transient pressure analysis of multiple-fractured horizontal wells in boxed reservoirs. Petroleum Exploration and Development 41(1): 74–78.

52.

Warren

Root

(1963) The behavior of naturally fractured reservoirs. Society of Petroleum Engineers Journal 3(3): 245–255.

53.

Wayne MA (2008) Geology of Carbonate Reservoirs: The Identification, Description, and Characterization of Hydrocarbon Reservoirs in Carbonate Rocks. Hoboken, NJ: Wiley.

54.

Wilkinson

(1989) New results for pressure transient behavior of hydraulically fractured wells. Paper 18950 presented at the SPE joint rocky mountain regional/low permeability reservoirs symposium and exhibition, Denver, Colorado, 6–8 March.

55.

Wang

Xing

, et al. (2017) Pressure-transient analysis for a vertically fractured well at an arbitrary azimuth in a rectangular anisotropic reservoir. Journal of Petroleum Science and Engineering 159: 279–294.

56.

Zeng

Zhao

Liu

(2012) A new model for reservoirs with a discrete-fracture system. Journal of Canadian Petroleum Technology 51(2): 127–136.

Pressure transient analysis of a vertical well with multiple etched fractures in carbonate reservoirs

Abstract

Keywords

Introduction

Mathematical model and algorithm

Point source function

Flow model in reservoir

Flow model in fracture

Single fracture solution

Multiple fractures solution

Model validation

Results and discussion

Flow regimes

Effect of number of fracture

Effect of fracture distribution

Effect of fracture conductivity

Effect of fracture storativity coefficient

Effect of interporosity flow parameter

Effect of different geometric shape of matrix blocks

Conclusions

Footnotes

Declaration of conflicting interests

Funding

References