Rate Decline Behavior of Selectively Completed Horizontal Wells in Naturally Fractured Oil Reservoirs

Abstract

Selectively completed horizontal wells (SCHWs) can significantly reduce cost of completing wells and delay water breakthrough and prevent wellbore collapse in weak formations. Thus, SCHWs have been widely used in petroleum development industry. SCHWs can shorten the effective length of horizontal wells and thus have a vital effect on production. It is significant for SCHWs to study their rate decline and flux distribution in naturally fractured reservoirs. In this paper, by employing motion equation, state equation, and mass conservation equation, three-dimension seepage differential equation is established and corresponding analytical solution is obtained by Laplace transform and finite cosine Fourier transform. According to the relationship of constant production and wellbore pressure in Laplace domain, dimensionless rate solution is gotten under constant wellbore pressure in Laplace domain. Dimensionless pressure and pressure derivate curves and rate decline curves are drawn in log-log plot and seven flow regimes are identified by Stehfest numerical inversion. We compared the simplified results of this paper with the results calculated by Saphir for horizontal wells in naturally fractured reservoirs. The results showed excellent agreement. Some parameters, such as outer boundary radius, storativity ratio, cross-flow coefficient, number and length of open segments, can obviously affect the rate integral and rate integral derivative log-log curves of the SCHWs. The proposed model in this paper can help better understand the flow regime characteristics of the SCHWs and provide more accurate rate decline analysis of the SCHWs data to evaluate formation.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In the early age, since horizontal wells can increase productivity, they became a popular method to develop oil and gas. Compared with vertical wells, horizontal wells can control severe water or gas coning problems, increase the connecting area with the reservoir, and reduce wellbore turbulence [1].

Although horizontal wells show a number of advantages, increasing wellbore length may lead to production imbalance along the wellbore, which can lead to water coning and decreased production. Certainly, the uneven rate distribution can lead to bottom-water break through [2]. However, horizontal open holes may be completed by employing prepacked screens because of their low cost. At the same time, prepacked screens can effectively minimize sand production. With the production of oil or gas, the sand will accumulate around the horizontal wellbore, which makes fluids of formation flow into wellbore by a section of wellbores opened and causes productivity decline. In order to reduce cost of completing horizontal well, delay water breakthrough, and prevent wellbore collapse in weak formations, SCHWs were used widely and actual production also proves the effectiveness of the method. An important feature of SCHWs was that only some segments of the wellbore are open to the formation. It is also observed that even if the entire length of the horizontal well is open or perforated, only some segments produce fluid [3].

To analyze the wellbore pressure and rate response of SCHWs, some engineers tend to use an effective horizontal well length to replace the open length of the horizontal well. This treatment assumes that the open length of the horizontal well is continuous instead of interval distribution. An analytical model was developed in real domain to predict the inflow performance of SCHWs and selectively completed vertical wells (SCVWs) [4, 5]. Their model considers the distribution of the open intervals. Kamal et al. [6] presented an analytical model of SCHWs by employing instantaneous-point-source solution and the superposition principle. They used this model to analyze actual cases where pressure transient is available. Yildiz [7–9] and Seyide [10] presented a model of SCHWs and SCVWs and derived asymptotic approximations of the model in Laplace space. Pressure and pressure-derivate log-log curves were plotted, flow characteristics were discussed, and each segment rate of SCHWs was analyzed. A new semianalytical model for predicting the performance of horizontal wells which were completed by inflow control devices in bottom-water reservoirs was presented. The coupled solution is developed for predicting the performance of horizontal wells in a box-shaped reservoir with bottom-water drive [11].

In order to analyze the rate decline curve of SCHWs, a mathematical model considering difference between horizontal and vertical permeability of SCHWs is established in naturally fractured reservoirs. Based on point source and the superposition principle, pressure analytical solution of the SCHWs under the condition of constant production in impermeable top and bottom boundary and lateral impermeable boundary by Laplace transform and finite cosine Fourier transform. Log-log curves of pressure and pressure-derivate and rate decline are drawn in naturally fractured reservoirs by employing Stehfest numerical inversion. Seven flow regimes, according to the characteristic of pressure-derivate curve, are identified and every flow regime characteristic is described in detail. This paper discusses that relevant parameters (storativity ratio, flow coefficient, number and length completed horizontal sections, etc.) have effect on pressure and rate decline curves. Corresponding solutions can be useful in completion design and rate decline in field practice.

2. Physical Model of SCHWs and Assumption

Horizontal wells are located in naturally fractured reservoirs with impermeable top and bottom boundary and lateral impermeable boundary and parallel to the upper and lower impermeable boundary. Horizontal well consists of $N_{\mathrm{N}\mathrm{c}}$ open segments and $N_{\mathrm{c}}$ completed segments (shown in Figure 1). The naturally fractured reservoir is structured by matrix system and natural fracture system. The dual-porosity media are assumed as the Warren-Root model and pseudo-steady cross flow exits between the matrix and fractures. Fluids flow into natural fractures from matrix firstly and flow into horizontal wellbore from natural fractures secondly by open segments (shown in Figure 2). In order to make mathematical model more reliable and accurate, some assumptions are listed as follows:

(a)

The fluid flow in the reservoir obeys Darcy’s law and law of isothermal percolation.

(b)

Flow is single phase and the fluid has constant and small compressibility and constant viscosity.

(c)

Formation permeability is anisotropic with three major directional permeability $k_x$, $k_y$, $k_z$.

(d)

Formation initial pressure is $p_{\mathrm{e}}$ and horizontal produced at a constant surface flow rate $q_{\mathrm{s}\mathrm{c}}$.

(e)

Horizontal well consists of $N_{\mathrm{N}\mathrm{c}}$ open segments and $N_{\mathrm{c}}$ completed segments, and fluids flow into wellbore by only open segments.

(f)

The length of the open segments and completed segments may be unequal, and each open segment may have a different skin effect and production rate.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

In this paper, we follow the point source theory adopted by Gringarten and Ramey [12] and Ozkan and Raghavan [13, 14] in order to obtain wellbore pressure under constant-rate production and rate distribution under constant-pressure production in naturally fractured reservoirs (see Notations section).

3. Mathematical Model of SCHWs

3.1. Point Source Model

3.1.1. Governing Equation

(A) Fracture System. The 3D governing equation describing transient fluid flow in natural fracture system can be written as follows:$$\begin{array}{}\text{(1)}& \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{f}}}{\partial r^{\mathrm{2}}}+\frac{\mathrm{1}}{r}\frac{\partial p_{\mathrm{f}}}{\partial r}+\frac{k_{\mathrm{f}\mathrm{v}}}{k_{\mathrm{f}\mathrm{h}}}\frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{f}}}{\partial z^{\mathrm{2}}}=\frac{{\varphi }_{\mathrm{f}}{C}_{\mathrm{f}\mathrm{t}}\mu }{k_{\mathrm{f}\mathrm{h}}}\frac{\partial p_{\mathrm{f}}}{\partial t}+\frac{{\varphi }_{\mathrm{m}}{C}_{\mathrm{m}\mathrm{t}}\mu }{k_{\mathrm{f}\mathrm{h}}}\frac{\partial p_{\mathrm{m}}}{\partial t}\end{array}$$

where $r=\sqrt{(x-x_{\mathrm{w}}{)}^2+(y-y_{\mathrm{w}}{)}^2}$

(B) Matrix System. The 3D governing equation describing transient fluid flow in natural fracture system can be written as follows:$$\begin{array}{}\text{(2)}& \frac{\alpha k_{\mathrm{m}}}{k_{\mathrm{f}\mathrm{h}}}\left(p_{\mathrm{m}}-p_{\mathrm{f}}\right)+\frac{{\varphi }_{\mathrm{m}}{C}_{\mathrm{m}}\mu }{k_{\mathrm{f}\mathrm{h}}}\frac{\partial p_{\mathrm{m}}}{\partial t}=\mathrm{0}\end{array}$$

3.1.2. Initial Conditions

The initial pressure is assumed to be equal and is represented by original formation pressure in naturally fractured oil reservoirs; thus$$\begin{array}{}\text{(3)}& p_{\mathrm{f}}\left(r,z,t=0\right)=p_{\mathrm{e}}\text{(4)}& p_{\mathrm{m}}\left(r,z,t=0\right)=p_{\mathrm{e}}\end{array}$$

3.1.3. Inner and Outer Boundary Conditions

It is assumed that production rate of point source is $\tilde{q}(t)$; thus inner boundary condition can be written as$$\begin{array}{}\text{(5)}& \underset{\epsilon \to \mathrm{0}}{\lim }{{\displaystyle \int }}_{z_{\mathrm{w}}-\epsilon /2}^{z_{\mathrm{w}}+\epsilon /2}\left[\underset{\delta \to 0}{\lim }2\pi \frac{k_{\mathrm{f}\mathrm{h}}}{\mu }{\left(r\frac{\partial p_{\mathrm{f}}}{\partial r}\right)}_{r=\epsilon }\right]dz=\tilde{q}\left(t\right)\end{array}$$

Corresponding outer boundary conditions can be expressed as for a laterally impermeable boundary, top, and bottom boundaries being(6) ∂ p f r , z = 0 , t ∂ z = 0 (7) ∂ p f r , z = h , t ∂ z = 0 (8) ∂ p f r = R e , z , t ∂ r = 0

3.2. Dimensionless Point Source Model in Laplace Domain

According to dimensionless variables definition in Table 1, (1)–(8) can be transformed into $$\begin{array}{c}\text{(9)}& \frac{{\partial }^2{p}_{\mathrm{f}\mathrm{D}}}{\partial r_{\mathrm{D}}^2}+\frac{1}{r_{\mathrm{D}}^{}}\frac{\partial p_{\mathrm{f}\mathrm{D}}}{\partial r_{\mathrm{D}}^{}}+\frac{{\partial }^2{p}_{\mathrm{f}\mathrm{D}}}{\partial z_{\mathrm{D}}^2}=\left(\omega \frac{\partial p_{\mathrm{f}\mathrm{D}}}{\partial t_{\mathrm{D}}^{}}+\left(1-\omega \right)\frac{\partial p_{\mathrm{m}\mathrm{D}}}{\partial t_{\mathrm{D}}^{}}\right)\lambda \left(p_{\mathrm{m}\mathrm{D}}-p_{\mathrm{f}\mathrm{D}}\right)+\left(1-\omega \right)\frac{\partial p_{\mathrm{m}\mathrm{D}}}{\partial t_{\mathrm{D}}^{}}=0\\ p_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{D}}=0\right)=0\\ p_{\mathrm{m}\mathrm{D}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{D}}=0\right)=0\\ \frac{\partial p_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},z_{\mathrm{D}}=0,t_{\mathrm{D}}\right)}{\partial z_{\mathrm{D}}^{}}=0\\ \frac{\partial p_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},z_{\mathrm{D}}=h_{\mathrm{D}},t_{\mathrm{D}}\right)}{\partial z_{\mathrm{D}}^{}}=0\\ \frac{\partial p_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}}=R_{\mathrm{e}\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{D}}\right)}{\partial r_{\mathrm{D}}^{}}=0\\ \underset{{\epsilon }_{\mathrm{D}}\to 0}{\lim }{{\displaystyle \int }}_{z_{\mathrm{w}\mathrm{D}}-{\epsilon }_{\mathrm{D}}/2}^{z_{\mathrm{w}\mathrm{D}}+{\epsilon }_{\mathrm{D}}/2}\left[\underset{{\delta }_{\mathrm{D}}\to 0}{\lim }{\left(r_{\mathrm{D}}^{}\frac{\partial p_{\mathrm{f}\mathrm{D}}}{\partial r_{\mathrm{D}}^{}}\right)}_{r_{\mathrm{D}}={\epsilon }_{\mathrm{D}}}\right]d{z}_{\mathrm{D}}=-h_{\mathrm{D}}{\tilde{q}}_{\mathrm{D}}\left(t_{\mathrm{D}}\right)\end{array}$$

Table 1

Dimensionless variables definition.

Variables	Dimensionless definition
Dimensionless pressure of fracture system	$p_{\mathrm{f}\mathrm{D}}=\frac{\mathrm{2}\pi {\stackrel{\hspace{0.17em}}{k}}_{\mathrm{f}\mathrm{h}}h}{q_{\mathrm{s}\mathrm{c}}\mu }\left(p_e-p_{\mathrm{f}}\right)$
Dimensionless pressure of matrix system	$p_{\mathrm{m}\mathrm{D}}=\frac{\mathrm{2}\pi k_{\mathrm{f}\mathrm{h}}h}{q_{\mathrm{s}\mathrm{c}}\mu }\left(p_e-p_{\mathrm{m}}\right)$
Dimensionless wellbore pressure	$p_{\mathrm{w}\mathrm{D}}=\frac{\mathrm{2}\pi k_{\mathrm{f}\mathrm{h}}h}{q_{\mathrm{s}\mathrm{c}}\mu }\left(p_e-p_{\mathrm{w}}\right)$
Dimensionless production time	$t_{\mathrm{D}}=\frac{k_{\mathrm{f}\mathrm{h}}t}{\mu {\left(\varphi C_{\mathrm{t}}\right)}_{\mathrm{f}+\mathrm{m}}{L}_{\mathrm{r}\mathrm{e}\mathrm{f}}^{\mathrm{2}}}$
Dimensionless distance	$r_{\mathrm{D}}=\frac{r}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}$
Dimensionless reservoir thickness	$h_{\mathrm{D}}=\frac{h}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}$
Dimensionless radius of impermeable circle boundary	$R_{\mathrm{e}\mathrm{D}}=\frac{R_{\mathrm{e}}}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}$
Dimensionless coordinate	$x_{\mathrm{D}}=\frac{x}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}{y}_{\mathrm{D}}=\frac{y}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}{z}_{\mathrm{D}}=\frac{z}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}\sqrt{\frac{k_{\mathrm{f}\mathrm{h}}}{k_{\mathrm{f}\mathrm{v}}}}$
Dimensionless x-y-z coordinate of point source	$x_{\mathrm{w}\mathrm{D}}=\frac{x_{\mathrm{w}}}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}{y}_{\mathrm{w}\mathrm{D}}=\frac{y_{\mathrm{w}}}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}{z}_{\mathrm{w}\mathrm{D}}=\frac{z_{\mathrm{w}}}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}\sqrt{\frac{k_{\mathrm{f}\mathrm{h}}}{k_{\mathrm{f}\mathrm{v}}}}$
Dimensionless length of open segment	$L_{\mathrm{N}\mathrm{c}\mathrm{D}}=\frac{L_{Nc}}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}$
Dimensionless wellbore radius	$r_{\mathrm{w}\mathrm{D}}=\frac{r_{\mathrm{w}}}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}$
Dimensionless mid-point of ith open segment	$x_{\mathrm{m}\mathrm{D}}=\frac{x_{\mathrm{m}}}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}$
Dimensionless continuous production	${\tilde{q}}_{\mathrm{D}}=\frac{\tilde{q}}{q_{\mathrm{s}\mathrm{c}}}$
Dimensionless infinitesimal vertical distance	${\epsilon }_{\mathrm{D}}=\frac{\epsilon }{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}$
Dimensionless infinitesimal radial distance	${\delta }_{\mathrm{D}}=\frac{\delta }{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}$

For convenience in derivation, by adopting Laplace transform with respect to $t_{\mathrm{D}}$ and solving pressure of matrix, respectively, through (9), thus (9) can be expressed as follows in Laplace domain.$$\begin{array}{c}\text{(10)}& \frac{{\partial }^{\mathrm{2}}{\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}}{\partial r_{\mathrm{D}}^{\mathrm{2}}}+\frac{\mathrm{1}}{r_{\mathrm{D}}^{}}\frac{\partial {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}}{\partial r_{\mathrm{D}}^{}}+\frac{{\partial }^{\mathrm{2}}{\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}}{\partial z_{\mathrm{D}}^{\mathrm{2}}}=u{\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}\frac{\partial {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},z_{\mathrm{D}}=\mathrm{0},s\right)}{\partial z_{\mathrm{D}}^{}}=\mathrm{0}\\ \frac{\partial {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},z_{\mathrm{D}}=h_{\mathrm{D}},s\right)}{\partial z_{\mathrm{D}}^{}}=\mathrm{0}\\ \frac{\partial {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}}=R_{\mathrm{e}\mathrm{D}},z_{\mathrm{D}},s\right)}{\partial r_{\mathrm{D}}^{}}=\mathrm{0}\\ \underset{{\epsilon }_D\to \mathrm{0}}{\lim }{{\displaystyle \int }}_{z_{\mathrm{w}\mathrm{D}}-{\epsilon }_{\mathrm{D}}/\mathrm{2}}^{z_{\mathrm{w}\mathrm{D}}+{\epsilon }_{\mathrm{D}}/\mathrm{2}}\left[\underset{{\delta }_{\mathrm{D}}\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}{\left(r_{\mathrm{D}}^{}\frac{\partial {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}}{\partial r_{\mathrm{D}}^{}}\right)}_{r_{\mathrm{D}}={\epsilon }_{\mathrm{D}}}\right]d{z}_{\mathrm{D}}=-h_{\mathrm{D}}{\stackrel{-}{\stackrel{~}{q}}}_{\mathrm{D}}\left(s\right)\end{array}$$

where $u=\left(\left(\lambda +s\omega (1-\omega )\right)/\left(\lambda +s(1-\omega )\right)\right)s$

By finite cosine transform with respect to $z_{\mathrm{D}}$, $(10)$ can be written as$$\begin{array}{c}\text{(11)}& \frac{d^{\mathrm{2}}{\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}}{d{r}_{\mathrm{D}}^{\mathrm{2}}}+\frac{\mathrm{1}}{r_{\mathrm{D}}^{}}\frac{d\hspace{0.17em}{\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}}{d{r}_{\mathrm{D}}^{}}=\left(u+\frac{n^2{\pi }^2}{h_{\mathrm{D}}^{\mathrm{2}}}\right){\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}\frac{d\hspace{0.17em}{\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}}=R_{\mathrm{e}\mathrm{D}},n,s\right)}{\partial r_{\mathrm{D}}^{}}=\mathrm{0}\\ \underset{{\delta }_D\to \mathrm{0}}{\lim }{\left(r_{\mathrm{D}}^{}\frac{d\hspace{0.17em}{\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}}{d{r}_{\mathrm{D}}^{}}\right)}_{r_D={\epsilon }_D}=-{\stackrel{-}{\stackrel{~}{q}}}_D\left(s\right)h_{\mathrm{D}}\cos \left(\frac{n\pi z_{\mathrm{w}\mathrm{D}}}{h_{\mathrm{D}}}\right)\end{array}$$

In deriving $(10)$, we have used the following finite cosine transform and inverse finite cosine transform.$$\begin{array}{}\text{(12)}& {\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},n,s\right)={{\displaystyle \int }}_{\mathrm{0}}^{h_{\mathrm{D}}}{\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},z_{\mathrm{D}},s\right)\cos \left(\frac{n\pi z_{\mathrm{D}}}{h_{\mathrm{D}}}\right)d{z}_{\mathrm{D}}\text{(13)}& {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},z_{\mathrm{D}},s\right)=\frac{\mathrm{1}}{h_{\mathrm{D}}}{\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},\mathrm{0},s\right)+\frac{\mathrm{2}}{h_{\mathrm{D}}}{\displaystyle \sum }_{n=\mathrm{1}}^{\infty }{\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}\left(r_{\mathrm{D}},n,s\right)\cos \left(\frac{n\pi z_{\mathrm{D}}}{h_{\mathrm{D}}}\right)\end{array}$$

3.3. Model Solution of Point Source

Equation $(10)$ is zero-order Bessel equation and inner and outer boundary conditions; general solution of zero-order Bessel equation can be written as$$\begin{array}{}\text{(14)}& {\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}=A{K}_{\mathrm{0}}\left({\epsilon }_{\mathrm{n}}{r}_{\mathrm{D}}\right)+B{I}_{\mathrm{0}}\left({\epsilon }_{\mathrm{n}}{r}_{\mathrm{D}}\right)\end{array}$$

where ${\epsilon }_{\mathrm{n}}=\sqrt{u+n^{\mathrm{2}}{\pi }^{\mathrm{2}}/h_{\mathrm{D}}^{\mathrm{2}}}$ (n=0,1,2,3…).

According to the properties of modified Bessel’s functions and outer boundary condition, the coefficient B can be expressed by$$\begin{array}{}\text{(15)}& B=\frac{K_{\mathrm{1}}\left({\epsilon }_n{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I_1\left({\epsilon }_n{R}_{\mathrm{e}\mathrm{D}}^{}\right)}A\end{array}$$

Hence, (14) also can be written as$$\begin{array}{}\text{(16)}& {\stackrel{⌢}{\stackrel{-}{p}}}_{\mathrm{f}\mathrm{D}}=A\left[K_{\mathrm{0}}\left({\epsilon }_{\mathrm{n}}{r}_{\mathrm{D}}\right)+\frac{K_1\left({\epsilon }_n{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I_{\mathrm{1}}\left({\epsilon }_n{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I}_0\left({\epsilon }_{\mathrm{n}}{r}_{\mathrm{D}}\right)\right]\end{array}$$

Combining with inner boundary condition, the coefficient A in (16) can be determined as follows.$$\begin{array}{}\text{(17)}& A={\stackrel{-}{\stackrel{~}{q}}}_{\mathrm{D}}{h}_{\mathrm{D}}\cos \left(\frac{n\pi z_{\mathrm{w}\mathrm{D}}}{h_{\mathrm{D}}}\right)\end{array}$$

Substituting (17) into (16) and employing finite cosine inverse transform, finally, (16) can be written as follows:$$\begin{array}{}\text{(18)}& {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}={\stackrel{-}{\stackrel{~}{q}}}_{\mathrm{D}}\left[K_{\mathrm{0}}\left({\epsilon }_{\mathrm{0}}{r}_{\mathrm{D}}\right)+\frac{K_{\mathrm{1}}\left({\epsilon }_{\mathrm{0}}{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I_{\mathrm{1}}\left({\epsilon }_{\mathrm{0}}{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I}_{\mathrm{0}}\left({\epsilon }_{\mathrm{0}}{r}_{\mathrm{D}}\right)+\mathrm{2}{\displaystyle \sum }_{n=\mathrm{1}}^{\infty }\left(K_{\mathrm{0}}\left({\epsilon }_{\mathrm{n}}{r}_{\mathrm{D}}\right)+\frac{K_{\mathrm{1}}\left({\epsilon }_n{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I_{\mathrm{1}}\left({\epsilon }_n{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I}_{\mathrm{0}}\left({\epsilon }_{\mathrm{n}}{r}_{\mathrm{D}}\right)\right)\cos \left(\frac{n\pi z_{\mathrm{D}}}{h_{\mathrm{D}}}\right)\cos \left(\frac{n\pi z_{\mathrm{w}\mathrm{D}}}{h_{\mathrm{D}}}\right)\right]\end{array}$$

Equation (18) is the point source solution in naturally fractured oil reservoirs.

3.4. Model Solution of Line Source for SCHWs

Taking the SCHWs shown in Figure 1, for example, although rate of the horizontal well in different location, open segments can be seen as a uniform rate horizontal line source for SCHWs. Thus we take the ith open segment as our research object, taking mid-point of ith open segment as origin coordinate. So we can get line source of $i$th open segment by integrating with respect to $x_{\mathrm{D}}$ from $x_{\mathrm{m}\mathrm{D}\mathrm{i}}$-$L_{\mathrm{N}\mathrm{c}\mathrm{D}\mathrm{i}}$/2 to $x_{\mathrm{m}\mathrm{D}\mathrm{i}}$-$L_{\mathrm{N}\mathrm{c}\mathrm{D}\mathrm{i}}$/2 for point source in (18). Finally$$\begin{array}{}\text{(19)}& {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}i}=\frac{{\stackrel{-}{q}}_{\mathrm{D}i}}{L_{\mathrm{N}\mathrm{c}\mathrm{D}i}}{{\displaystyle \int }}_{x_{\mathrm{m}\mathrm{D}i}-L_{\mathrm{N}\mathrm{c}\mathrm{D}i}/\mathrm{2}}^{x_{\mathrm{m}\mathrm{D}i}+L_{\mathrm{N}\mathrm{c}\mathrm{D}i}/\mathrm{2}}\left[K_{\mathrm{0}}\left({\epsilon }_{\mathrm{0}}\sqrt{{\left(x_{\mathrm{D}}-\alpha \right)}^{\mathrm{2}}}\right)+\frac{K_{\mathrm{1}}\left({\epsilon }_{\mathrm{0}}{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I_{\mathrm{1}}\left({\epsilon }_{\mathrm{0}}{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I}_{\mathrm{0}}\left({\epsilon }_{\mathrm{0}}\sqrt{{\left(x_{\mathrm{D}}-\alpha \right)}^{\mathrm{2}}}\right)+\mathrm{2}{\displaystyle \sum }_{n=\mathrm{1}}^{\infty }\left(K_{\mathrm{0}}\left({\epsilon }_{\mathrm{n}}\sqrt{{\left(x_{\mathrm{D}}-\alpha \right)}^{\mathrm{2}}}\right)+\frac{K_{\mathrm{1}}\left({\epsilon }_n{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I_{\mathrm{1}}\left({\epsilon }_n{R}_{\mathrm{e}\mathrm{D}}^{}\right)}{I}_{\mathrm{0}}\left({\epsilon }_{\mathrm{n}}\sqrt{{\left(x_{\mathrm{D}}-\alpha \right)}^{\mathrm{2}}}\right)\right)\cos \left(\frac{n\pi z_{\mathrm{D}}}{h_{\mathrm{D}}}\right)\cos \left(\frac{n\pi z_{\mathrm{w}\mathrm{D}}}{h_{\mathrm{D}}}\right)\right]\end{array}$$

It is noted that (19) is only valid to compute the pressure anywhere. We use ${\mathrm{z}}_{\mathrm{D}}$=${\mathrm{z}}_{\mathrm{w}\mathrm{D}}$+${\mathrm{r}}_{\mathrm{w}\mathrm{D}}$ to calculate wellbore surface pressure. Dimensionless wellbore radius for an anisotropic reservoir is given as follows.$$\begin{array}{}\text{(20)}& r_{\mathrm{w}\mathrm{D}}=\frac{r_{\mathrm{w}}}{L_{\mathrm{r}\mathrm{e}\mathrm{f}}}\sqrt{{\left(\frac{k}{k_z}\right)}^{\mathrm{0.25}}+{\left(\frac{k}{k_z}\right)}^{-\mathrm{0.25}}}\end{array}$$

According to geometric relations shown in Figure 3, mid-point of ith open segment can be determined as follows.$$\begin{array}{}\text{(21)}& x_{\mathrm{m}\mathrm{D}i}=L_{\mathrm{c}\mathrm{D}i}+L_{\mathrm{N}\mathrm{c}\mathrm{D}i}+\frac{L_{\mathrm{c}\mathrm{D}\mathrm{1}}}{\mathrm{2}}\end{array}$$

[figure omitted; refer to PDF]

With (19) and by applying the superposition principle, the pressure responses at point ($x_{\mathrm{D}}$, $y_{\mathrm{D}}$, $z_{\mathrm{D}}$) caused by all segments can be obtained as follows:$$\begin{array}{}\text{(22)}& {\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}}\left(x_{\mathrm{D}},\mathrm{0},z_{\mathrm{D}}\right)={\displaystyle \sum }_{i=1}^N{\stackrel{-}{p}}_{\mathrm{f}\mathrm{D}i}\left({\beta }_{\mathrm{D}},{\beta }_{\mathrm{D}i}\right)={\displaystyle \sum }_{i=1}^N{\stackrel{-}{q}}_{\mathrm{D}i}{F}_{\mathrm{D}i}\left({\beta }_{\mathrm{D}},{\beta }_{\mathrm{D}i}\right)\end{array}$$

where ${\beta }_{\mathrm{D}}=\left(x_{\mathrm{D}},\mathrm{0},z_{\mathrm{D}}\right)$; ${\beta }_{\mathrm{D}i}=\left(x_{\mathrm{D}i},\mathrm{0},z_{\mathrm{D}i}\right)$

However, it is also required that the sum of the flow rates for each open segment be equal to the total flow rate; that is,$$\begin{array}{}\text{(23)}& {\displaystyle \sum }_{i=1}^N{\stackrel{-}{q}}_{\mathrm{D}i}=\frac{\mathrm{1}}{s}\end{array}$$

Combining with (19) and (22), then the matrix form can be formulated as follows.$$\begin{array}{}\text{(24)}& \left[\begin{array}{ccccc}{\overline{F}}_{\mathrm{D}\mathrm{1}}\left({\beta }_{\mathrm{D}\mathrm{1}},{\beta }_{\mathrm{D}\mathrm{1}}\right)& {\overline{F}}_{\mathrm{D}\mathrm{2}}\left({\beta }_{\mathrm{D}\mathrm{1}},{\beta }_{\mathrm{D}\mathrm{2}}\right)& \dots & {\overline{F}}_{\mathrm{D}N}\left({\beta }_{\mathrm{D}\mathrm{1}},{\beta }_{\mathrm{D}N}\right)& -\mathrm{1}\\ {\overline{F}}_{\mathrm{D}\mathrm{1}}\left({\beta }_{\mathrm{D}\mathrm{2}},{\beta }_{\mathrm{D}\mathrm{1}}\right)& {\overline{F}}_{\mathrm{D}\mathrm{2}}\left({\beta }_{\mathrm{D}\mathrm{2}},{\beta }_{\mathrm{D}\mathrm{2}}\right)& \dots & {\overline{F}}_{\mathrm{D}N}\left({\beta }_{\mathrm{D}\mathrm{2}},{\beta }_{\mathrm{D}N}\right)& -\mathrm{1}\\ \dots & \dots & \dots & \dots & \dots \\ {\overline{F}}_{\mathrm{D}\mathrm{1}}\left({\beta }_{\mathrm{D}\mathrm{M}},{\beta }_{\mathrm{D}\mathrm{1}}\right)& {\overline{F}}_{\mathrm{D}\mathrm{2}}\left({\beta }_{\mathrm{D}\mathrm{2}},{\beta }_{\mathrm{D}\mathrm{M}}\right)& \dots & {\overline{F}}_{\mathrm{D}N}\left({\beta }_{\mathrm{D}N},{\beta }_{\mathrm{D}N}\right)& -\mathrm{1}\\ \mathrm{1}& \mathrm{1}& \dots & \mathrm{1}& \mathrm{0}\end{array}\right]\times \left[\begin{array}{c}{\overline{q}}_{\mathrm{D}\mathrm{1}}\\ {\overline{q}}_{\mathrm{D}\mathrm{2}}\\ \dots \\ {\overline{q}}_{\mathrm{D}N}\\ {\stackrel{-}{p}}_{\mathrm{w}\mathrm{D}}\end{array}\right]=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \dots \\ \mathrm{0}\\ \frac{\mathrm{1}}{s}\end{array}\right]\end{array}$$

The dimensionless wellbore flow rate for the constant-pressure production in naturally fractured reservoirs can be determined by dimensionless pressure with the constant-rate production in the Laplace domain [15]:$$\begin{array}{}\text{(25)}& {\stackrel{-}{q}}_{\mathrm{D}}=\frac{\mathrm{1}}{s^{\mathrm{2}}{\stackrel{-}{p}}_{\mathrm{w}\mathrm{D}}}\end{array}$$

In order to be consistent with current literature, we use the Fetkovich [16] definitions of the dimensionless decline variables ($t_{\mathrm{D}\mathrm{d}}$ and $q_{\mathrm{D}\mathrm{d}}$) which are given below. The $t_{\mathrm{D}\mathrm{d}}$ function is given in terms of dimensionless variables as$$\begin{array}{}\text{(26)}& t_{\mathrm{D}\mathrm{d}}=\frac{\mathrm{2}}{R_{\mathrm{e}\mathrm{D}}^{\mathrm{2}}}\frac{t_{\mathrm{D}}}{\ln R_{\mathrm{e}\mathrm{D}}^{}-\mathrm{0.5}}\end{array}$$

In a similar fashion, the $q_{\mathrm{D}\mathrm{d}}$ function is given in terms of dimensionless variables as$$\begin{array}{}\text{(27)}& q_{\mathrm{D}\mathrm{d}}=\left[\ln R_{\mathrm{e}\mathrm{D}}^{}-\mathrm{0.5}\right]q_{\mathrm{D}}\end{array}$$

The rate integral and rate integral derivative functions introduced by McCray [17] are given in dimensionless form below. The dimensionless rate integral function, $q_{\mathrm{D}\mathrm{i}}$, is given as$$\begin{array}{}\text{(28)}& q_{\mathrm{D}\mathrm{d}\mathrm{i}}=\frac{N_{\mathrm{p}\mathrm{D}\mathrm{d}}}{t_{\mathrm{D}\mathrm{d}}}=\frac{\mathrm{1}}{t_{\mathrm{D}\mathrm{d}}}{{\displaystyle \int }}_0^{t_{\mathrm{D}\mathrm{d}}}{q}_{\mathrm{D}\mathrm{d}}\left(x\right)dx\end{array}$$

And the dimensionless rate integral derivative function, $q_{\mathrm{D}\mathrm{d}\mathrm{i}\mathrm{d}}$, is given as$$\begin{array}{}\text{(29)}& q_{\mathrm{D}\mathrm{d}\mathrm{i}\mathrm{d}}=q_{\mathrm{D}\mathrm{d}\mathrm{i}\mathrm{d}}-q_{\mathrm{D}\mathrm{d}}\end{array}$$

4. Model Verification

To verify the model and solutions derived in the above section, a relatively particular case is considered and pressure and pressure-derivate curves generated by our solution are compared to well-test stimulator Saphir. Fluid flow into wellbore is treated as infinite conductivity, but rate distribution in wellbore is no-uniform. Therefore, based on different dimensionless variable definition between this paper and well-test stimulator, we can set ω=0.2 for this paper and well-test stimulator, λ=0.01 for this paper, λ=6.25 × 10⁻⁸ with $L_{\mathrm{r}\mathrm{e}\mathrm{f}}$=40 for well-test stimulator. Other parameters can be set as ω=0.2, $L_{\mathrm{h}}$=400, $L_{\mathrm{N}\mathrm{c}}$=40, $L_{\mathrm{c}}$=0, $N_{\mathrm{c}}$=0, $N_{\mathrm{N}\mathrm{c}}$=10. The comparisons presented in Figure 4 suggest that the results calculated by our model are consistent with that obtained by the well-test simulator, which verifies the credibility of the model presented in this article. At the same time, by analyzing each segment rate distribution in different time, we can confirm imbalanced fluid inflow along the wellbore (see Figure 5).

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

5. Discussions and Analysis

5.1. Flow Regimes

In order to study the flow regimes of SCHWs in naturally fractured oil reservoirs more graphically, type curves of pressure response and production rate performance are illustrated in Figures 6 and 7 by Stehfest [18] numerical inversion. According to the dimensionless pressure-derivate characteristic, pressure response curves of SCHWs in naturally fractured oil reservoirs are divided into seven flow regimes and the important basic data is shown in Table 2.

Table 2

Important basic data for SCHWs.

Parameters (unit)	Value
Wellbore radius (m)	0.1
Outer boundary radius (m)	10000
Horizontal well length (m)	400
Length of each open segment (m)	25
Length of each completed segment (m)	100
Number of open segments (dimensionless)	4
Number of completed segments (dimensionless)	3
Reference length (m)	40
Storativity ratio (dimensionless)	0.2
Flow coefficient (dimensionless)	0.01

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

Period I is the first radial (FR) flow period. During this period, the flow regime is radial flow around open segment in vertical direction (see Figure 8(a)). The pressure derivative curve is a horizontal line with a value of “1/(4$L_{\mathrm{h}\mathrm{D}}$).” Curves of rate, rate integral, and integral derivative exhibit a downward line. Because each length of open segment is different, first radial flow may not occur. If open segment is very small, spherical flow can appear during this period.

[figures omitted; refer to PDF]

Period II is first linear (FL) flow period, in which fluid flow in the reservoir is parallel to the upper and lower boundary of the reservoir and each open segment is independent during first linear flow stage (see Figure 8(b)). The main characteristic of pressure-derivative is a line with a half slope in this stage. Corresponding rate integral derivative curve becomes gentle in this stage. Similarly, when open segment is shorter compared with complete segment, first linear flow cannot be appearing.

Period III is second pseudo-radial (SPR) flow, in which the pressure derivative curve is horizontal line of “0.5/$N_{\mathrm{N}\mathrm{c}}$.” This flow period is exhibited when the lateral distance between open segments is relatively large. Before and during this flow period, each open segment has its own drainage area and behaves independently without interference from open segment (see Figure 8(c)).

Period IV is second linear (SL) flow. Pressure wave propagates to drainage area controlled by each open segment, and interference between open segments occurs. After the superposition of the pressure waves, pressure waves propagate continually as time goes. The second linear flow can be formed in natural fracture reservoir (see Figure 8(d)). The pressure derivative curve during this period is exhibited as a one-half slope straight line again and characteristics of rate integral derivative curve do not appear in this stage.

Period V is the cross-flow stage, in which fluid flows into natural fracture from matrix firstly when the SCHWs are put into production. The pressure of natural fracture system will gradually decrease, causing pressure difference between natural fracture system and matrix system. Because existence of pressure drop between natural fracture system and matrix system led to cross-flow from natural fracture system to matrix system, the characteristic of pressure derivative during stage is “dip.” Corresponding rate integral derivative curve also exhibits a “dip” in this stage.

Period VI is late pseudo-radial (LPR) flow stage. After cross-flow flow stage, the pressures in natural fracture system and matrix system gradually incline to equilibrium. Pseudo-radial flow around SCHWs is formed in naturally fractured reservoirs (see Figure 8(e)). Pressure derivative exhibits a horizontal line of “0.5” during pseudo-radial flow in log-log plot. Corresponding rate integral derivative curve is also a slanted line.

Period VII is characteristic of closed boundary. Pressure waves propagate to circular impermeable outer boundary during this stage. Curves of pressure derivative exhibit unite-slope line and corresponding rate integral and derivative curve coincide and exhibit negative unite-slope line.

5.2. Sensitivity of Parameters

Figure 9 shows the effect of outer boundary radius on dimensionless rate integral $q_{\mathrm{D}\mathrm{d}\mathrm{i}}$ and rate integral derivative $q_{\mathrm{D}\mathrm{d}\mathrm{i}\mathrm{d}}$. We can know that outer boundary radius has effect on whole flow regime. With increase of outer boundary radius, the value of dimensionless rate integral and rate integral derivative curves is smaller, which indicates that larger outer boundary radius can lead to the smaller rate decline curves in whole flow regime.

[figure omitted; refer to PDF]

Figure 10 shows the effect of storativity ratio on dimensionless rate integral $q_{\mathrm{D}\mathrm{d}\mathrm{i}}$ and rate integral derivative $q_{\mathrm{D}\mathrm{d}\mathrm{i}\mathrm{d}}$. It is obvious that storativity ratio mainly has significant effect on FR, FL, SPR, and cross-flow regime. The smaller storativity ratio represents that storativity ability of naturally fractured reservoirs is smaller. It can be clearly observed that the smaller the value of storativity ratio, the deeper and wider the “dip” in rate integral derivative curve during cross-flow. In addition, smaller value of storativity ratio leads to smaller rate integral of FR, FL, and SPR.

[figure omitted; refer to PDF]

Figure 11 shows the effect of cross-flow coefficient on dimensionless rate integral $q_{\mathrm{D}\mathrm{d}\mathrm{i}}$ and rate integral derivative $q_{\mathrm{D}\mathrm{d}\mathrm{i}\mathrm{d}}$. It is obvious that flow coefficient mainly has significant effect on SPR, cross-flow, and LPR regime. The larger cross-flow coefficient represents that cross-flow ability from natural fracture to matrix is larger. It can be clearly observed that the smaller the value of cross-flow coefficient is, the later the “dip” in rate integral derivative curve during cross-flow appears. Duration of LPR regime is shorter.

[figure omitted; refer to PDF]

Figures 12 and 13 show the effect of number of open (or completed) segments on dimensionless rate integral $q_{\mathrm{D}\mathrm{d}\mathrm{i}}$ and rate integral derivative $q_{\mathrm{D}\mathrm{d}\mathrm{i}\mathrm{d}}$. Though number of open (or completed) segments is different, it is assumed that length of SCHWs is equal. Number of completed segments has an effect on SPR, cross-flow, and LR regimes mainly. With increase of completed segment, rate of fluid flow into wellbore decreases under constant wellbore pressure. Therefore, higher number of completed segments can lead to larger value of rate integral derivative curves in log-log plot (see Figure 12). At the same time, increasing number of completed segments can delay water breakthrough and prevent wellbore collapse; it can also lead to small rate for every open segment, which makes total rate decrease under constant wellbore pressure (see Figure 13).

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

Figures 14 and 15 show the effect of length of completed segment on dimensionless rate integral $q_{\mathrm{D}\mathrm{d}\mathrm{i}}$ and rate integral derivative $q_{\mathrm{D}\mathrm{d}\mathrm{i}\mathrm{d}}$. It is assumed that lengths of SCHWs and numbers of completed segments are equal, while completed segment length is different. It is obvious that length of completed segment has an effect on each flow regime. Compared with LR regime, length of completed segment has an obvious influence on rate integral derivative curve during FR, FL, SPR, and cross-flow regime. Longer length of completed segment leads to low rate integral derivative, which is caused by larger pressure loss from formation to wellbore (see Figure 14). With the increase of length of completed segment, open segment becomes more and more short. As the time of production continues to increase, longer completed segment makes each open segment rate smaller (see Figure 15).

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

6. Conclusion

In this work, we have developed a solution to compute the rate decline of SCHWs with constant wellbore pressure. According to characteristic of pressure-derivative curves under constant production and rate integral derivative curve under constant wellbore pressure, flow regimes of SCHWs are identified. Specific conclusions can be drawn as follows.

$(1)$ An analytical model is proposed in this paper to obtain rate decline response and flux distribution of SCHWs in naturally fractured reservoirs under constant wellbore pressure. Pressure transient responses and Blasingame rate decline curve are generated and discussed.

$(2)$ In addition, the seven flow periods observed for SCHWs mainly include first radial flow, first linear flow, second pseudo-radial flow, second pseudo-radial flow, second linear flow, late pseudo-radial flow stage, and characteristic of closed boundary.

$(3)$ By comparing results of simplified model in this paper with the results calculated by Saphir for horizontal well in naturally fractured reservoirs, the results showed excellent agreement. Imbalanced flux distribution along the wellbore is verified again.

$(4)$ The model illustrated how the dimensionless rate integral and rate integral derivative log-log curves are influenced by some parameters (such as outer boundary radius, storativity ratio, cross-flow coefficient, and number and length of open segments).

$(5)$ The proposed model in this paper can be used to interpret rate decline signals more accurately for SCHWs in naturally fractured oil reservoirs and provide more accurate dynamic parameters which are important for efficient reservoir development.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

You-jie Xu and Qi-guo Liu contributed equally to this work (co-first authors).

Acknowledgments

This article was supported by the National Major Research Programme for Science and Technology of China (Grant No. 2017ZX05009-004 and No. 2016ZX05015-003).

Glossary

Notations

$C_{\mathrm{f}\mathrm{t}}$ :Total compressibility of natural fracture system and oil, atm⁻¹

$C_{\mathrm{m}\mathrm{t}}$ :Total compressibility of matrix system and oil, atm⁻¹

$h$ :Reservoir thickness, cm

$k$ : Equivalent permeability, μm², $k=\sqrt[\mathrm{3}]{k_x{k}_y{k}_z}$

$k_{\mathrm{f}\mathrm{h}}$ :Horizontal permeability of natural fracture system, μm²

$k_{\mathrm{f}\mathrm{v}}$ :Vertical permeability of natural fracture system, μm²

$k_{\mathrm{m}}$ :Permeability of matrix system, μm²

$k_{\mathrm{z}}$ :Vertical permeability, μm²

$L_{\mathrm{c}i}$ :Length of $i$th completed segment, cm

$L_{\mathrm{N}\mathrm{c}i}$ :Length of $i$th open segment, cm

$L_{\mathrm{r}\mathrm{e}\mathrm{f}}$ :Reference length, cm

$N$ :Open segment number, dimensionless

$N_{\mathrm{p}}$ :Cumulative production, cm³

$p_{\mathrm{e}}$ :Initial reservoirs pressure, atm

$p_{\mathrm{f}}$ :Pressure of natural fracture system, atm

$p_{\mathrm{m}}$ :Pressure of natural matrix system, atm

$p_{\mathrm{w}}$ :Wellbore pressure of natural matrix system, atm

$q$ :Production under constant wellbore pressure, cm³

$q_{\mathrm{d}}$ :Decline rate function as defined by Fetkovich, cm³

$q_{\mathrm{d}\mathrm{i}}$ :Decline rate integral as defined by McCray, cm³

$q_{\mathrm{d}\mathrm{i}\mathrm{d}}$ :Decline rate integral derivative function as defined by McCray

$q_{\mathrm{s}\mathrm{c}}$ :Production rate under the standard conditions, cm³/s

$\tilde{q}(t)$ :Surface production rate of a point source, cm³/s

$r$ :Radial distance, cm

$R_{\mathrm{e}}$ :Radius of impermeable circle boundary, cm

$r_{\mathrm{w}}$ :Wellbore radius, cm

$s$ :Laplace variables

$t$ : $[\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\hspace{0.17em}\hspace{0.17em}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}$ , s

$t_{\mathrm{d}}$ :Decline time, s

$x$ :x-coordinates, cm

$x_{\mathrm{m}}$ :Mid-point of $i$th open segment, cm

$x_{\mathrm{w}}$ :x-coordinates of a point source, cm

$y$ :y-coordinates, cm

$y_{\mathrm{w}}$ :y-coordinates of a point source, cm

$z$ :z-coordinates, cm

$z_{\mathrm{w}}$ :z-coordinates of a point source, cm

$\alpha$ :Shape factor of dual-porosity system, cm⁻²

$\delta$ :Infinitesimal radial distance, cm

$\epsilon$ :Infinitesimal vertical distance, cm

$\lambda$ :Cross-flow coefficient of dual-porosity reservoirs, dimensionless

$\mu$ :Viscosity at current reservoir pressure, cp

$\phi$ :Reservoir porosity, dimensionless

$\omega$ :Storativity ratio of dual-porosity reservoirs, dimensionless

$I_{\mathrm{0}}(x)$ :The first kind modified Bessel function, zero order

$K_{\mathrm{0}}(x)$ :The second kind modified Bessel function, zero order

$I_{\mathrm{1}}(x)$ :The first kind modified Bessel function, first order

$K_{\mathrm{2}}(x)$ :The second kind modified Bessel function, first order.

$\mathrm{D}$ :Dimensionless

$\mathrm{f}$ :Natural fracture system

$\mathrm{m}$ : Matrix system.

_:Laplace domain

$\stackrel{⌢}{\hspace{0.17em}}$ :Finite cosine transform.

References

[1] R. Leon-Ventura, G. Gonzalez-G., H. Leyva-G., "Evaluation of Horizontal Well Production," Proceedings of the SPE International Petroleum Conference and Exhibition in Mexico,DOI: 10.2118/59062-MS, .

[2] M. M. Saggaf, "A vision for future upstream technologies," Journal of Petroleum Technology, vol. 60 no. 3, pp. 54-98, 2008.

[3] F. Brons, V. Marting, "The Effect of Restricted Fluid Entry on Well Productivity," Journal of Petroleum Technology, vol. 13 no. 02, pp. 172-174, DOI: 10.2118/1322-G-PA, 2013.

[4] P. A. Goode, D. J. Wilkinson, "Inflow performance of partially open horizontal wells," Journal of Petroleum Technology, vol. 43 no. 8, pp. 983-987, DOI: 10.2118/19341-PA, 1991.

[5] L. Larsen, "The Pressure-Transient Behavior of Vertical Wells With Multiple Flow Entries," Proceedings of the SPE Annual Technical Conference and Exhibition,DOI: 10.2118/26480-MS, .

[6] M. Kamal, I. Buhidma, S. Smith, W. Jones, "Pressure-transient analysis for a well with multiple horizontal sections," Proceedings of the SPE Annual Technical Conference and Exhibition,DOI: 10.2118/26444-MS, .

[7] T. Yildiz, E. Ozkan, "Transient pressure behaviour of selectively completed horizontal wells," Proceedings of the SPE Annual Technical Conference and Exhibition,DOI: 10.2118/28388-MS, .

[8] T. Yildiz, Y. Cinar, "Inflow Performance and Transient Pressure Behavior of Selectively Completed Vertical Wells," SPE Reservoir Engineering, vol. 1 no. 5, pp. 467-473, 1998.

[9] T. Yildiz, "Productivity of horizontal wells completed with screens," SPE Reservoir Evaluation and Engineering, vol. 7 no. 5, pp. 342-350, DOI: 10.2118/88997-PA, 2004.

[10] H. O. Seyide, "Performance Analysis of a Selectively Completed Horizontal Well," Proceedings of the SPE Nigeria Annual International Conference and Exhibition,DOI: 10.2118/189057-MS, .

[11] W. Luo, H.-T. Li, Y.-Q. Wang, J.-C. Wang, "A new semi-analytical model for predicting the performance of horizontal wells completed by inflow control devices in bottom-water reservoirs," Journal of Natural Gas Science and Engineering, vol. 27, pp. 1328-1339, DOI: 10.1016/j.jngse.2015.03.001, 2015.

[12] A. C. Gringarten, H. J. Ramey, "The use of source and green’s function in solving unsteady-flow problem in reservoir," SPE Journal, vol. 13 no. 5, pp. 285-296, DOI: 10.2118/3818-PA, 1973.

[13] E. Ozkan, R. Raghavan, "New solutions for well-test-analysis problems. Part 1. Analytical considerations," SPE Formation Evaluation, vol. 6 no. 3, pp. 359-368, DOI: 10.2118/18615-PA, 1991.

[14] E. Ozkan, R. Raghavan, "New solutions for well-test-analysis problems. Part 2. Computational considerations and applications," SPE Formation Evaluation, vol. 6 no. 3, pp. 369-378, DOI: 10.2118/18616-PA, 1991.

[15] A. F. van Everdingen, W. Hurst, "The application of the Laplace transformation to flow problems in reservoirs," Journal of Petroleum Technology, vol. 1 no. 12, pp. 305-324, DOI: 10.2118/949305-G, 1949.

[16] M. Fetkovich, "Decline curve analysis using type curves," Journal of Petroleum Technology, vol. 32 no. 6, pp. 1065-1077, DOI: 10.2118/4629-pa, 2013.

[17] T. L. Mccray, Reservoir analysis using production decline data and adjusted time, 1990.

[18] H. Stehfest, "Numerical inversion of Laplace transform," Communications of the ACM, vol. 13 no. 1, pp. 47-49, DOI: 10.1145/361953.361969, 1970.