Introduction

A fractured medium is formed by the effect of stresses that break the rock that contains tectonic fractures between blocks of rock, and there is no fluid interchange between rock and fractures. These are type I reservoirs in accordance with Nelson (2001) classification of naturally fractured reservoirs.

Some examples are fractured igneous rock or fractured reservoirs classified as type I by Nelson (2001), with extension fractures where displacement is perpendicular to the walls of the fracture. When these fractures are filled with hydrocarbons, they are called fissures. In the field of fracture mechanics, it is common to classify them, as mode I fractures. Figure 1 shows a mode fracture with a perpendicular displacement to its walls.

Fig. 1
figure 1

Mode I fracture (Fossen 2010)

Full size image

There are various types of extension fractures, such as fissures, joints, and veins, which are observed in outcrops, cores, and fractured reservoirs.

Figure 2 displays different outcrops with calcite-filled extension fractures in limestone and joints in sandstone. These parallel fractures present several apertures and connected fracturing patterns.

Fig. 2
figure 2

Extension fractures. a Calcite-filled extension fractures in limestone. b Joints in sandstone. Image Courtesy United States Geological Survey; image source: Earth Science World Image Bank (AGI). Marli Miller, University of Oregon. http://www.earthscienceworld.org/images

Full size image

Figure 3 shows a limestone core with an open tectonic fracture. These types of fractures, considered as extension fractures, present a visible perpendicular displacement. Some models of hydrocarbon reservoirs have used these kinds of natural fractures.

Fig. 3
figure 3

Tectonic fractures in a core of NFTR

Full size image

The theory of fluid flow in fractured media developed by Barenblatt and Kochina (1960) is based on the assumption of constant rock properties. Barenblatt’s model consisted of two media, matrix and fractures, with matrix–fracture transfer which would generate a pressure gradient during hydrocarbon production. Despite the medium being fractured, fluid flow is represented using Darcy’s law.

The proposed analytical techniques assume constant rock properties, which yield a constant diffusivity. In the present paper, we use a Navier–Stokes solution called the Couette’s equation to model fluid flow in extension fractures.

In practice, fractured reservoirs types I, II, and III classified by Nelson (2001) are non-stress sensitive media. The analysis has been developed chiefly with the aim of obtaining analytical expressions for the solution of the mathematical flow model, for naturally fractured tectonic reservoirs.

An analytical model for non-stress-sensitive naturally fractured tectonic reservoirs is developed; it is solved analytically using Cole–Hopf transform for the case of an infinite reservoir and Couette’s flow that includes a quadratic gradient term.

The study of the present paper is based on outcrops, core samples, and field data; we show that NFTRs could be modeled using Couette’s flow, considering the effects of the nonlinear gradient term to describe fluid flow.

Darcy and Couette equation

Darcy’s law is frequently used and sometimes unknowing its basic assumptions. The most restrictive application condition is related to the Reynolds number; namely, that fluid flow is dominated by viscous forces, considering laminar flow for Reynolds number, Re, which means a number smaller than unity (Muskat 1946).

Various authors give different limiting values for Darcy’s laminar flow, between a range of Re from 3 to 10 (Polubarinova-Kochina 1962). However, Muskat (1946) discussed that Darcy’s law can be applied to reservoirs flow problems whose conditions yield Reynolds number smaller that unity.

Figure 4 shows different Reynolds numbers for applicability of Darcy’s law, considering linear nonlinear laminar flow and turbulent flow. We applied the Couette’s flow in the nonlinear laminar zone with Reynolds numbers ranging between 5 and 13 (Couland et al. 1986).

Fig. 4
figure 4

Applicability of Darcy’s law (Virtual Campus in Hydrology and Water Resources Management 2014)

Full size image

The Reynolds number and the basic Darcy equation may be stated as:

$$Re = \frac{{\overline{{D_{\text{p}} }} v\rho }}{\mu }$$
(1)

For natural fractures, it can be expressed as:

$$Re = \frac{qa\rho }{\mu A\phi }$$
(1a)

where

$$v = - \frac{k\rho }{\mu }\nabla\Phi$$
(2)

where \(\Phi , p/\rho + gz\); \(\Phi ,{\text{flow}}\,{\text{potential}}\); \(p,\,{\text{formation}}\,{\text{pressure}}\); \(\rho ,\,{\text{oil}}\, {\text{density}}\); \(g,\,{\text{gravity}}\); \(z,\,{\text{elevation}}\); \(\mu ,{\text{oil}}\,{\text{dynamic}}\,{\text{viscosity}}\); A, 2πrh; \(A,\,{\text{area}}\); \(r,\,{\text{radius}}\); \(h,{\text{thickness}}\,{\text{formation}}\); \(k,{\text{total}}\,{\text{permeability}}\); \(Re,\,{\text{Reynolds}}\,{\text{number}}\); \(\overline{{D_{\text{p}} }} ,\,{\text{average}}\,{\text{pore}} \,{\text{diameter}}\); \(a,\,{\text{fracture}}\,{\text{aperture}}\); \(v_{\text{d}} ,\,{\text{specific}}\,{\text{discharge}}\); \(\phi ,\,{\text{porosity}}\).

Darcy’s law is valid for the median of the flow probability distribution (Scheidegger 1960) and is based on the assumption that fluid flow is inertialess.

It can be stated that for a heterogeneous, anisotropic and fractured porous medium, the upper limit critical Reynolds number for laminar flow ranges from 0.1 to 10. The transition to quadratic flow (without reaching turbulence), see the nonlinear laminar section of Fig. 4, was demonstrated by Schneebeli (1955). The nonlinear seepage flow law will be parabolic at Re > 13, with deviation from linearity (Barenblatt et al. 1990; Couland et al. 1986); also nonlinear corrections to Darcy’s law at low Reynolds numbers for periodic porous media have been described (Firdaouss et al. 1997). Values of Reynolds number between 5 and 13 were calculated with numerical experiments based on the Navier–Stokes equations (Couland et al. 1986; Stark 1972).

Nonlinear flow is found in fractured porous media. Consequently, the Couette equation can be used for analytical modeling because it has a quadratic flow profile that is an exact solution for the Navier–Stokes equation; this equation is similar to cubic law and/or Boussinesq’s formula. The cubic law estimates the fluid flow rate for flow through fractures systems; usually, this equation is used in naturally fractured tectonic reservoirs (NFTRs), considering the laminar flow of a viscous fluid between parallel flat plates (Barros-Galvis et al. 2015; Potter and Wiggert 2007). On the other hand, Singh and Sharma (2001) used an extension of the three dimensional Couette flow to study the channel flow and the effect of the permeability of the porous medium.

The application of Couette or Darcy equations is associated with the Reynolds number. Figure 5 shows the high-velocity fluid flow for naturally fractured tectonic reservoirs, which is related to the Reynolds number, too. For radial flow in a reservoir, two zones will be observed, a high-velocity zone of radius, r hv, and other low velocity zone for greater radii that r hv.

Fig. 5
figure 5

Stabilized zone of non-Darcy flow for radial flow toward a well (high velocity)

Full size image

For radial flow, it has been described that: for a flowing well the high-velocity flow stabilizes at a radius, which the Reynolds number is one. Namely, linear laminar flow and Darcy flow are reached.

The red circle represents the inner (minimum) radius for Darcy’s flow; for \(r < r_{\text{hv}}\), flow is under high-velocity conditions, and Couette equation is used, which Reynolds number is greater than unity.

We can derive the seepage law, using the Navier–Stokes equations by means of integration (Barenblatt et al. 1990) and Couette equation. In this paper, we use and discuss Couette equation to describe fluid flow in natural fractures.

Analytical model

In order to develop this mathematical model, some considerations are as follows:

  1. 1.

    Fluid is stored and transported in natural tectonic fractures.

  2. 2.

    Single phase. Flow of an undersaturated oil reservoir, so that the fluid is a liquid (Craft and Hawkins 1991).

  3. 3.

    Porosity, permeability, and density rock are constants. So, they do not depend on either stresses or fluid pressure.

  4. 4.

    Isotropic permeability

  5. 5.

    Liquid is uncompressible or slightly compressible, in consequence fluid density changes exponentially with respect to pressure (Muskat 1946).

  6. 6.

    Isothermal fluid flow of small and constant compressibility.

Analytical modeling

The analytical model is based on a partial differential equation that describes the fluid flow in fractures. In developing this equation, we combine: a continuity equation or law of conservation of mass, a flow law such as the Couette’s equation, and an equation of state (Barros-Galvis 2015).

A linear diffusivity equation depicting the flow of incompressible liquid in a fractured medium can be obtained.

For homogeneous media, the flow law used is Darcy’s law (Pedrosa 1986; Chin and Raghavan 2000; Marshall 2009), satisfying the Reynolds number Re < 1.

Figure 6 shows a tectonic fracture represented as two parallel surfaces. The flow between these plates is taken to be in the x direction, and since there is no flow in the y direction, pressure will only be a function of the x direction. In addition, there are no inertia, viscous, or external forces in the y direction.

Fig. 6
figure 6

Flow profile between two inclined parallel surfaces (Potter  and Wiggert 2007)

Full size image

Fluid flow in an extension fracture is modeled using an exact solution to the Navier–Stokes equation, referred to as the general Couette flow (see Eq. 3). This equation describes fluid flow through extension tectonic fractures (Currie 2003):

$$u\left( y \right) = - \frac{1}{2\mu }\frac{{{\text{d}}\left( {p_{\text{f}} + \gamma H} \right)}}{{{\text{d}}x}}y\left( {a - y} \right) + \frac{U}{a}y$$
(3)

where \(p_{\text{f}} ,\,{\text{fracture}}\,{\text{pressure}}\); \(u(y),\,{\text{velocity}}\,{\text{profile}}\); \(U,\,{\text{upper }}\,{\text{surface}}\,{\text{velocity}}\); \(\beta ,\,{\text{inclination}}\,{\text{angle}}\); \(y,\,{\text{vertical}}\,{\text{direction}}\); \(x,\,{\text{horizontal}}\, {\text{direction}}\); \(\gamma ,\,{\text{specific}}\,{\text{weight}}\); \(H,\,{\text{vertical}}\, {\text{distance}}\); \(a,\,{\text{fracture}}\, {\text{aperture}}\); \(\mu ,\,{\text{oil}}\, {\text{dynamic}}\,{\text{viscosity}}\).

Equation (3) shows that the velocity profile across the flow field is parabolic. There are two ways of inducing flow between two parallel surfaces: (1) applying a pressure gradient and (2) the upper surface moves in the x direction with constant velocity U.

In this paper, we induce flow applying a pressure gradient; the maximum velocity occurs in y = a/2, so that the application of a pressure gradient presupposes that the upper surface will be fixed, and flow can be described by Poiseuille equation; in consequence, the Poiseuille flow is a specific case of the general Couette flow. Saidi (1987) used the Poiseuille equation for the flow in channel and fractures.

The use of the maximum velocity at y = a/2 (Fig. 6) in Couette’s equation indicates the highest fluid flow rate into a discontinuity, for a fixed pressure gradient. In the solution that follows, gravity is neglected, and Eq. (3) can be written as:

$$u\left( y \right) = - \frac{1}{2\mu }\frac{{{\text{d}}p_{\text{f}} }}{{{\text{d}}x}}y\left( {a - y} \right)$$
(4)

Equation (4) corresponds to the steady flow pressure distribution through two inclined parallel surfaces; it can be used to derive an equation for the flow rate using the following expression:

$$q = \int u{\text{d}}A$$
(4a)

Substituting Eq. (4) into Eq. (4a):

$$\begin{aligned} q & = \mathop \int \limits_{0}^{a} \frac{1}{2\mu }\left( {y^{2} - ay} \right)\frac{{\partial \left( {p_{\text{f}} + \gamma H} \right)}}{\partial x}{\text{d}}y \\ & = - \frac{{a^{3} }}{12\mu }\frac{{\partial \left( {p_{\text{f}} + \gamma H} \right)}}{\partial x} \\ \end{aligned}$$
(4b)

Average velocity, \(\bar{u}(y)\), is defined as \(\bar{u}(y) = q /A\), where cross-sectional area A = a * L, L = 1 (unitary length), and q is oil flow rate.

$$\bar{u}\left( y \right) = - \frac{{a^{2} }}{12\mu }\frac{{\partial \left( {p_{\text{f}} + \gamma H} \right)}}{\partial x}$$
(4c)

Equation (4b) is known as the cubic law, and in accordance with the flow direction, the sign may be positive or negative. If fractures are horizontal, then H = 0.

Equation (4) is similar to Darcy’s equation; considering y = a/2, this equation can be rewritten, for the maximum velocity profile, u(y)max:

$$u\left( y \right)_{\hbox{max} } = - \frac{{a^{2} }}{8\mu }\frac{{{\text{d}}p_{\text{f}} }}{{{\text{d}}x}} = - \frac{{a^{2} }}{8\mu }\nabla p_{\text{f}}$$
(5)

Fracture permeability k f can be expressed as follows (Aguilera 1995):

$$k_{\text{f}} = 8.35 \times 10^{6 } \,(a^{2} ) \,{\text{darcys}}$$
(6)

In Eq. (6), the aperture (a) is in centimeters. These equations combine Poiseuille’s law for capillary flow and Darcy’s law for flow of liquids in permeable beds. Craft and Hawkins (1991) used Eq. (6) to estimate permeability in channels or smooth fractures surfaces, with constant wide aperture.

Singha and Al-Shammeli (2012) reported and estimated values for tectonic and non-tectonic fractures conductivity that were matched with well testing for field cases. In consequence, fractures are not smooth surfaces. They used Poiseuille’s law, calibrating hydraulic conductivity, and it is given by:

$$C_{\text{f}} = k_{\text{f}} h = \frac{{a^{3} }}{{12 \times 0.98 \times 10^{ - 6 } }}$$
(7)

where C f is tectonic fracture conductivity in md m and a is fracture aperture in mm:

$$a = \sqrt[3]{{ k_{\text{f}} \times 11.76 \times 10^{ - 6 } }}$$
(7a)

Equation (7a) is developed for diffuse fractures and fracture corridors, with aperture values ranging from 0.0618744 and 39.9288 mm, respectively.

Equations (7) and (7a) are determined for rough and realistic limestone fractures in a tight carbonate reservoir in the Middle East with tectonic fractures, with fracture porosity of 0.2%. This reservoir is classified as type I, with low fracture porosity and high permeability.

Equation (6) applies to smooth and ideal fractures. Normally, reservoir modeling assumes ideal fractures, generating uncertainty. Substituting Eq. (6) into Eq. (5):

$$u\left( y \right) = - \frac{{k_{\text{f}} }}{{66.8 \times 10^{6 } \mu }}\nabla p_{\text{f}}$$
(8)

Equation (8) describes fluid velocities for smooth and ideal fractures, where a is fracture aperture in cm and \(k_{\text{f}}\) is fracture permeability in Darcy’s.

Applying Eq. (7), the calculated conductivity values is 20 md m for 6.24 × 10−2 mm average fracture aperture in 1 m of formation thickness, h. Equation (7a) can be substituted into Eq. (5):

$$u\left( y \right) = - \frac{{\left( { k_{\text{f}} \times 11.76 \times 10^{ - 6 } } \right)^{{{\raise0.7ex\hbox{$2$} \!\mathord{\left/ {\vphantom {2 3}}\right.\kern-0pt} \!\lower0.7ex\hbox{$3$}}}} }}{8\mu } \nabla p_{\text{f}}$$
$$u\left( y \right) = - \frac{{\left( { k_{\text{f}} } \right)^{{{\raise0.7ex\hbox{$2$} \!\mathord{\left/ {\vphantom {2 3}}\right.\kern-0pt} \!\lower0.7ex\hbox{$3$}}}} }}{{32.26 \times 10^{8 } \mu }} \nabla p_{\text{f}}$$
(8a)

It can be observed two constants, \(C = 66.8 \times 10^{6 }\) for Eq. (8), and C = 32.26 × 108 for Eq. (8a).

Equations (8) and (8a) are expressions similar to Darcy’s law, where u(y) = v:

$$u\left( y \right) = - \frac{{k_{\text{f}} }}{C\mu }\nabla p_{\text{f}}$$
(9)

To derive a partial differential equation for fluid flow in a fractured medium, we should combine a flow law with the continuity equation (Matthews and Russell 1967; Lee et al. 2003).

The continuity equation can be expressed using a derivative or integral equation, which they are equivalent; considering the former case, this equation is given by Eq. (10):

$$\frac{{\partial \left( {\rho \phi_{\text{f}} } \right)}}{\partial t} = - \nabla \cdot \left( {\rho v} \right)$$
(10)

where \(v,\,{\text{Couette}}\text{'}{\text{s}}\, {\text{velocity}}\,{\text{equation}}\); \(\phi ,\,{\text{fracture}}\,{\text{porosity}}\); \(\rho ,\,{\text{oil}} \,{\text{density}}\); \(t,\,{\text{time}}\).

Substituting Eq. (9) into Eq. (10) gives:

$$\frac{{\partial \left( {\rho \phi_{\text{f}} } \right)}}{\partial t} = - \nabla \cdot \left[ {\rho \left( { - \frac{{k_{\text{f}} }}{C\mu }\nabla p_{\text{f}} } \right)} \right]$$
(11)

Equation (11) can be expressed as

$$\frac{{\partial \left( {\rho \phi_{\text{f}} } \right)}}{\partial t} = \nabla \left( {\frac{{\rho k_{\text{f}} }}{{\mu^{*} }}} \right)\nabla p_{\text{f}} + \left( {\frac{{\rho k_{\text{f}} }}{{\mu^{*} }}} \right)\nabla^{2} p_{\text{f}}$$
(12)

where

$$\mu^{*} = C\mu$$

Each of the two terms on the right-hand side of Eq. (12) involves the permeability, viscosity, and porosity, which are constants. However, fluid density is pressure dependent. The present problem is restricted to single-phase liquids and slightly compressible liquids with constant compressibility, c, defined by Eq. (13):

$$c = \frac{1}{\rho }\frac{{{\text{d}}\rho }}{{{\text{d}}p_{\text{f}} }}$$
(13)

The liquid (oil) compressibility is a dominant term in total system compressibility (see Eq. 22);

For constant compressibility c, integration of Eq. (13) gives

$$\rho = \rho_{\text{o}} {\text{e}}^{{c\left( {p_{\text{f}} - p_{\text{i}} } \right)}}$$
(14)

where \(\rho\) is oil density; and ρ o is considered at initial pressure, p i, and p f is a reference pressure in the fracture.

The derivative of Eq. (14) with respect to pressure yields Eq. (15):

$$\frac{\partial \rho }{{\partial p_{\text{f}} }} = \left[ {\rho_{\text{o}} {\text{e}}^{{{\text{c}}\left( {p_{\text{f}} - p_{\text{i}} } \right)}} } \right]c = \rho c$$
(15)

Applying the chain rule and substituting Eq. (15):

$$\frac{\partial \rho }{\partial t} = \frac{\partial \rho }{\partial t}\frac{{\partial p_{\text{f}} }}{{\partial p_{\text{f}} }} = \frac{\partial \rho }{{\partial p_{\text{f}} }}\frac{{\partial p_{\text{f}} }}{\partial t} = \rho c\frac{{\partial p_{\text{f}} }}{\partial t}$$
(16)

For one-dimensional flow, the gradient in the right-hand side of first term of Eq. (12) can be written as:

$$\nabla \left( {\frac{{\rho k_{\text{f}} }}{{\mu^{*} }}} \right) = \frac{\partial }{\partial x}\left( {\frac{{\rho k_{\text{f}} }}{{\mu^{*} }}} \right)\frac{{\partial p_{\text{f}} }}{{\partial p_{\text{f}} }} = \frac{\partial \rho }{{\partial p_{\text{f}} }}\left( {\frac{{k_{\text{f}} }}{{\mu^{*} }}} \right)\frac{{\partial p_{\text{f}} }}{\partial x}$$
$$\nabla \left( {\frac{{\rho k_{\text{f}} }}{{\mu^{*} }}} \right) = \rho c\left( {\frac{{k_{\text{f}} }}{{\mu^{*} }}} \right)\frac{{\partial p_{\text{f}} }}{\partial x} = \rho c\left( {\frac{{k_{\text{f}} }}{{\mu^{*} }}} \right)\nabla p_{\text{f}}$$
(17)

Substituting Eqs. (16) and (17) into Eq. (12):

$$\rho \phi_{\text{f}} c\frac{{\partial p_{\text{f}} }}{\partial t} = \rho c\left( {\frac{{k_{\text{f}} }}{{\mu^{*} }}} \right)\left( {\nabla p_{\text{f}} } \right)^{2} + \left( {\frac{{\rho k_{\text{f}} }}{{\mu^{*} }}} \right)\nabla^{2} p_{\text{f }}$$
(18)

Transposing terms, Eq. (18) may be expressed:

$$\frac{{\partial p_{\text{f}} }}{\partial t}\frac{1}{D} = \left[ {c\left( {\nabla p_{\text{f}} } \right)^{2} + \nabla^{2} p_{\text{f}} } \right]$$
(19)

where the hydraulic diffusivity D is defined by Eq. (20)

$$D = \frac{{k_{\text{f}} }}{{\phi_{\text{f}} c\mu^{*} }}$$
(20)

A similar procedure using Eq. (8a) gives, for rough fractures (rf)

$$D_{\text{rf}} = \frac{{k_{\text{f}}^{2/3} }}{{\phi_{\text{f}} c\mu^{*} }}$$
(20a)

Tectonic reservoirs with extension fractures present low matrix porosity and permeability. Permeability and effective porosity of fractures are dominant variables in fluid flow; consequently, formation porosity and permeability are approximated to fractures properties. In this paper, matrix properties are considered and included in total permeability and porosity of the reservoir. In other words, petrophysical properties of the reservoir are considered in fractures system.

Initial and boundary conditions in radial coordinate are:

  1. 1.

    p f = p i at t = 0 for all r.

  2. 2.

    \(\left( {r\partial p_{\text{f}} /\partial r} \right)_{{r_{\text{w}} }} = - 6q\mu /\pi kh\) for t > 0.

To develop the solution, this boundary condition is replaced by the line source condition:

$$\mathop {\lim }\limits_{r \to 0} \left( {r\partial p_{\text{f}} /\partial r} \right)_{{r_{\text{w}} }} = - 6q\mu /\pi kh \,\,{\text{for}}\,t > 0$$
  1. 3.

    p f(r, t) = p i as \(r \to \infty\) for all t.

Equation (19) is a nonlinear partial differential equation and can be referred to as a nonlinear diffusivity equation (see Eqs. 19, 20, and 20a). This equation represents an analytical model for non-stress-sensitive naturally fractured tectonic reservoir, which describes fluid flow in the fracture system for an oil fractured reservoir, considering a nonlinear term of quadratic gradient (∇p f)2, and without matrix–fracture transfer.

Many papers have been published for the single-phase flow in homogeneous reservoirs that do not include the nonlinear pressure gradient term in the diffusivity equation, considered as small negligible pressure gradient, constant rock properties, and a fluid of small and constant compressibility; in effect, the nonlinear quadratic term is usually neglected (see the first right-hand side term of Eq. 20) for liquid flow (the fluid compressibility c has a small value) (Samaniego et al. 1979; Dake 1998; Matthews and Russell 1967). In addition, for an infinite reservoir the wellbore pressure predicted by this linear Darcy solution may be significantly smaller than that predicted by the Couette solution at large times. On the other hand, Jelmert and Vik (1996) and Odeh and Babu (1988) concluded that the consideration of the nonlinear quadratic term gives results significantly smaller in pressure prediction and recommended its use as the use pressure solution; although this result was also demonstrated by (Chakrabarty et al. 1993) for wellbore pressure prediction for a closed outer boundary, the authors stated that the linear pressure solution is unsatisfactory and should be applied with caution, stating that an infinite reservoir has a 5% error for large dimensionless times.

Others papers have presented solutions for the nonlinear transient flow model including a quadratic gradient term by using transformations (Friedel and Voigt 2009; Aadnoy and Finjord 1996; Chakrabarty et al. 1993; Cao et al. 2004); however, they assumed a homogeneous porous medium.

Mathematical model and solution for constant rate radial flow in an infinite naturally fractured reservoir

Our aim is to apply a mathematical transformation to reduce a nonlinear equation to linear equation diffusivity, for a naturally fractured system.

The differences between Darcy and Couette equations applied to the linear diffusivity equation have been described. Previous authors have not included the nonlinear pressure gradient term in the nonlinear diffusivity equation for fractures, or homogeneous systems. In both cases, fluid flow equations (Darcy and Couette equations) are used in this solution, considering parallel (slab) and single fractures geometry.

The diffusivity equation models mass and momentum transfer in the reservoir. The phenomenological description for fluid flow in NFTR is given by: (1) complex diffusion in tectonic fractures and (2) hydrodynamics as a result of a pressure gradient in well.

Complex diffusion contains various types of diffusion: molecular diffusion, surface diffusion, Knudsen diffusion, and convection due to gradient pressure. The real fractured system is heterogeneous and anisotropic, and their diffusion processes depend on fractures aperture or porous diameter (Cunningham and Williams 1980; Treybal 1980). In consequence, fast complex diffusion is reached in a nonlinear laminar flow, which may be modeled by Couette equation. Finally, momentum transfer is modeled using Couette equation.

Hydrodynamics in the wellbore is governed by pressure gradient caused by fluid flow. The fluid velocity is related to oil production rate, and Couette or Darcy equation application depends on the value of the Reynolds number. When Reynolds number is greater that unity Couette equation is applied. This application impacts the boundary condition of the diffusivity equation.

For bulk and slab block fractures properties, the following expressions should be considered, assuming a quasi-impermeable or aphanitic matrix.

$$\phi = \phi_{\text{m}} + \phi_{\text{f}}$$
(21)
$$c = \frac{{c_{\text{o}} + c_{\text{w}} \phi_{\text{m}} + c_{\text{m}} \phi_{\text{m}} + c_{\text{f}} \phi_{\text{f}} }}{{\phi_{\text{f}} }}$$
(22)
$$k = \frac{{k_{\text{f}} \left[ {N\pi \left( {\frac{a}{2}} \right)^{2} } \right] + k_{\text{m}} \left[ {A - N\pi \left( {\frac{a}{2}} \right)^{2} } \right]}}{A}$$
(23)
$$k_{\text{slab}} = \frac{{k_{\text{f}} a }}{d}$$
(24)

where \(\phi ,\,{\text{total}}\,{\text{porosity}}\); \(\phi_{\text{m}} ,\,{\text{matrix}}\, {\text{porosity}}\); \(\phi_{\text{f}} ,\,{\text{fracture}}\, {\text{porosity}}\); \(k,\,{\text{total}}\, {\text{permeability}}\); \(k_{\text{m}} ,\,{\text{matrix}}\,{\text{permeability}}\); \(k_{\text{f}} ,\,{\text{fracture}}\,{\text{permeability}}\); \(c,\,{\text{total}} \,{\text{compresibility}}\); \(a,\,{\text{fracture}}\,{\text{aperture}}\); \(c_{\text{o}} ,\,{\text{oil}} \,{\text{compresibility}}\); \(c_{\text{m}} ,\,{\text{matrix}}\,{\text{compressibility}}\); \(c_{\text{w}} ,\,{\text{water}}\,{\text{compressibility}}\); \(c_{\text{f}} ,\,{\text{fracture}}\, {\text{compresibility}}\); \(d,\,{\text{distance}}\, {\text{between}}\,{\text{fractures}}\); \(N,\,{\text{number}}\, {\text{of}}\,{\text{fractures}}\,{\text{per}}\,{\text{flow }}\,{\text{section}}\); \(k_{\text{slab}} ,\,{\text{parallel}}\,{\text{fractures}}\,{\text{permeability}}\).

Equations (21), (22), (23), and (24) are used by Reiss (1980) and Aguilera (1995).

Case 1A: Mathematical linear model and its solution for constant rate radial flow in an infinite reservoir using Darcy equation

The proposed mathematical model, used and documented by van Everdingen and Hurst (1949) and Carslaw and Jaeger (1959), Matthews and Russell (1967), Earlougher (1977), Streltsova (1988), Dake (1998), and Lee et al. (2003), will be solved for the pressure behavior of a well producing under constant rate in a radial infinite reservoir, using the documented solution for the diffusivity equation.

The diffusivity equation for the radial flow in an infinite homogeneous reservoir, with its initial and boundary conditions (Matthews and Russell 1967; Lee et al. 2003), is given by Eq. (25):

$$\frac{{\partial^{2} p}}{{\partial r^{2} }} + \frac{1}{r}\frac{\partial p}{\partial r} = \frac{\phi \mu c}{k}\frac{\partial p}{\partial t}$$
(25)

Initial and boundary conditions:

  1. 1.

    p(r, 0) = p i at t = 0 for all r

  2. 2.

    \(\left( {r\partial p/\partial r} \right)_{{r_{\text{w}} }} = - q\mu /2\pi kh\quad{\text{for}}\,t > 0 .\)

To develop the solution, this boundary condition is replaced by the line source condition:

$$\mathop {\lim }\limits_{r \to 0} \left( {r\partial p/\partial r} \right)_{{r_{\text{w}} }} = - q\mu /2\pi kh\,\,\,\,{\text{for }}\,t > 0$$
  1. 3.

    p(r, t) = p i as \(r \to \infty\) for all t.

The solution for the reservoir pressure is given by:

$$p_{\rm i} - p\left( {r,t} \right) = \frac{q\mu }{4\pi kh}\left[ { - E_{\rm i} \left( { - \frac{{\phi \mu cr^{2} }}{4kt}} \right)} \right]$$
(26)

where E i is the exponential integral function.

$$- E_{\rm i} \left( { - x} \right) = \mathop \int \limits_{x}^{\infty } \frac{{{\text{e}}^{ - u} }}{u}{\text{d}}u$$
(27)

For x < 0.0025,

$$-E_{\rm i} \left({- x} \right) \cong - \ln \left(\varepsilon x \right) = \ln \left({1/x} \right) - 0.5772.$$

The \(\varepsilon\) symbol is Euler’s constant, equal to 1.78. Thus, for (4kt/ϕμcr 2) > 100

$$p\left( {r,t} \right) = p_{\rm i} - \frac{q\mu }{4\pi kh}\left[ {\ln \left( {\frac{kt}{{\phi \mu cr^{2} }}} \right) + 0.80907} \right]$$
(28)

Equations (26) and (28) are well known as the line source solution. The flowing wellbore pressure is expressed by Eq. (29):

$$p_{\text{wf}} = p_{\rm i} - \frac{q\mu }{4\pi kh}\left[ {\ln \left( {\frac{kt}{{\phi \mu cr_{\text{w}}^{2} }}} \right) + 0.80907} \right]$$
(29)

The solution presented in Eq. (29) describes the actual finite–wellbore infinite reservoir, based on the assumption of a small wellbore radius, where \(p,\,{\text{formation}}\,{\text{pressure}}\); \(p_{\text{wf}} ,\,{\text{wellbore}}\, {\text{pressure}}\); \(p_{\text{i}} ,\,{\text{initial}}\,{\text{formation}}\,{\text{pressure}}\); \(\phi ,\,{\text{total }}\,{\text{porosity}}\); \(r, \,radius\); \(r_{\text{w}} ,\,{\text{wellbore}}\,{\text{radius}}\); \(t,\,{\text{time}}\); \(k,\,{\text{total}}\,{\text{permeability}}\); \(c,\,{\text{total}}\, {\text{compressibility}}\); \(h,\,{\text{formation}}\, {\text{thickness}}\); \(\mu ,\,{\text{oil}}\, {\text{viscosity}}\); \(q,\,{\text{oil}}\,{\text{flow rate}}\).

Commonly, Darcy’s law is used to model the inner boundary condition for the corresponding radial flow model.

Case 1B

Linear mathematical model and its solution for constant rate radial flow in an infinite reservoir, using Couette flow in the diffusivity equation and in the inner boundary condition.

For nonlinear laminar or Couette flow, a similar flow equation to Eq. (25) is obtained, but the boundary condition is different. The internal boundary condition is the cubic law that is implicit in Couette’s flow equation which can be considered as a solution for the Navier–Stokes equations (Witherspoon 1980).

$$\frac{{\partial^{2} p_{\text{f}} }}{{\partial r^{2} }} + \frac{1}{r}\frac{{\partial p_{\text{f}} }}{\partial r} = \frac{{\phi \mu^{*} c}}{k}\frac{{\partial p_{\text{f}} }}{\partial t}$$
(30)

The initial and boundary conditions are:

  1. 1.

    p f(r, 0) = p i at t = 0 for all r

  2. 2.

    \(\left( {r\partial p_{\text{f}} /\partial r} \right)_{{r_{\text{w}} }} = - 6q\mu /\pi ha^{2} \,\,\,\,\,{\text{for}}\,t > 0.\)

To develop the solution, this boundary condition is replaced by the following condition (which is similar to the line source approximation for radial flow):

$$\begin{aligned} & \mathop {\lim }\limits_{r \to 0} \left( {r\partial p_{\text{f}} /\partial r} \right)_{{r_{\text{w}} }} = - 6q\mu /\pi ha^{2} \,\,\,{\text{for }}\,t > 0 \\ & {\text{where}} \\ & q = Au\left( y \right) \\ & A = 2\pi rh, \\ & {\text{and}} \\ & u\left( y \right) = - \left( {a^{2} /12\mu } \right)\nabla p_{\text{f}} , \,{\text{then}} \\ & q = - \pi rh\left( {a^{2} /6\mu } \right)\nabla p_{\text{f}} \\ \end{aligned}$$

where

  1. 3.

    p f(r, t) = p i as \(r \to \infty\) for all t.

The solution of Eq. (30) is:

$$p_{\text{f}} \left( {r,t} \right) - p_{\text{i}} = - \frac{3q\mu }{{\pi ha^{2} }}\left[ {E_{\text{i}} \left( {\frac{{\phi \mu^{*} cr^{2} }}{4kt}} \right)} \right]$$
(31)

where E i represents the exponential integral function; for argument values <0.0025 [(4kt/ϕμ*cr 2) > 100], the logarithmic approximation for the wellbore pressure is:

$$p_{\text{wf}} = p_{\text{i}} - \frac{3q\mu }{{\pi ha^{2} }}\left[ {\ln \left( {\frac{kt}{{\phi \mu^{*} cr_{\text{w}}^{2} }}} \right) + 0.80907} \right]$$
(32)

The resultant equations are given in “Appendix 1,” which also outlines the solution procedure.

Case 1C

Linear mathematical model and its solution for constant rate using Darcy’s law in diffusivity equation and the Couette equation as inner boundary condition.

When Darcy’s flow equation is used to describe the flow in the reservoir, this physical phenomenon may be presented when the Reynolds number near and far the well is less than or equal to unity.

$$\frac{{\partial^{2} p_{\text{f}} }}{{\partial r^{2} }} + \frac{1}{r}\frac{{\partial p_{\text{f}} }}{\partial r} = \frac{\phi \mu c}{k}\frac{{\partial p_{\text{f}} }}{\partial t}$$
(33)

Again, initial and boundary conditions are:

  1. 1.

    p f = p i at t = 0 for all r.

  2. 2.

    \(\left( {r\partial p_{\text{f}} /\partial r} \right)_{{r_{\text{w}} }} = - 6q\mu /\pi ha^{2} \quad{\text{for}}\,t > 0.\)

  3. 3.

    \(p_{\text{f}} \to p_{\text{i}}\) as \(r \to \infty\) for all t.

Applying the last procedure of Case 1B, the solution to Eq. (33) for constant rate and radial flow in an infinite reservoir is given by:

$$p_{\text{i}} - p_{\text{f}} \left( {r,t} \right) = \frac{3q\mu }{{\pi ha^{2} }}\left[ { - E_{\text{i}} \left( { - \frac{{\phi \mu cr^{2} }}{4kt}} \right)} \right]$$
(34)

For pressure at the wellbore, r = r w, and (4kt/ϕμcr 2) > 100, Eq. (34) can be written:

$$p_{\text{wf}} = p_{\text{i}} - \frac{3q\mu }{{\pi ha^{2} }}\left[ {\ln \left( {\frac{kt}{{\phi \mu cr_{\text{w}}^{2} }}} \right) + 0.80907} \right]$$
(35)

Solution of the nonlinear non-stress-sensitive partial differential equation

The Cole–Hopf transformation was employed to obtain a solution to the Burger’s nonlinear partial differential equation (Burgers 1974; Ames 1972). Also, it has been used to solve a nonlinear diffusion problem for the flow of compressible liquids through homogeneous porous media (Marshall 2009). The Cole–Hopf transformation is a mathematical technique, through which a nonlinear partial differential equation may be reduced to linear partial differential equation.

Equation (20) describes the fluid flow in a fractured non-stress sensitive reservoir, with a nonlinear pressure term. These reservoirs are well known as type I, or single fracture model according to (Nelson 2001; Cinco-Ley 1996) classifications, respectively.

Case 2

Flow in extension fractures, without matrix–fracture fluid transfer and nonlinear diffusion. Oil is stored in the extension fractures and oil production flows through them to the wellbore. There is no matrix–fractures transfer because matrix porosity and permeability is very low, or tends to 0%, so that this matrix does not practically contain fluids.

It was observed that the transformation y = F(p f) applied in Burgers equation generated a linear partial equation of type ∂y/∂t = D∇2 y, and this concept was utilized to solve the nonlinear diffusivity equation (Ames 1972; Burgers 1974; Marshall 2009):

$$\frac{{\partial p_{\text{f}} }}{\partial t} = D\nabla^{2} p_{\text{f }} + D\left[ {\frac{{F^{{\prime \prime }} \left( {p_{\text{f}} } \right)}}{{F^{{\prime }} \left( {p_{\text{f}} } \right)}}} \right]\left( {\nabla p_{\text{f}} } \right)^{2}$$
(36)

with a quadratic gradient term. It can be observed that Eq. (20) is similar or equivalent to Eq. (36). To find a solution to Eq. (20), we would solve Eq. (36) for F. Then

$$y = F\left( {p_{\text{f}} } \right) = \frac{1}{c}\left[ {{\text{e}}^{{\left( {cp_{\text{f}} + a} \right)}} } \right] + b$$
(37)
$$F^{\prime } \left( {p_{\text{f}} } \right) = {\text{e}}^{{\left( {cp_{\text{f}} + a} \right)}}$$
(38)
$$F^{\prime \prime } \left( {p_{\text{f}} } \right) = c{\text{e}}^{{\left( {cp_{\text{f}} + a} \right)}}$$
(39)

where a and b are arbitrary integration constants generated due to the integration of F′(p f) and F″(p f).

Equation (37) is named the Cole–Hopf transformation. If a = b = 0 (Tong and Wang 2005), the transformation variable vanishes at some reference pressure. Then

$$y = \frac{1}{c}{\text{e}}^{{\left( {cp_{\text{f}} } \right)}} \Leftrightarrow p_{\text{f}} = \frac{1}{c} \ln (cy)$$
$$\frac{{\partial p_{\text{f}} }}{\partial y} = \frac{1}{cy}$$
(40)
$$\frac{{\partial^{2} p_{\text{f}} }}{{\partial y^{2} }} = - \frac{1}{{cy^{2} }}$$
(41)

Our goal is to eliminate \(\left( {\nabla p_{\text{f}} } \right)^{2}\); to accomplish this task, we need to derive expressions for: ∂p f/∂t, \(\nabla^{2} p_{\text{f }}\) and (∇p f)2:

$$\frac{{\partial p_{\text{f}} }}{\partial t} = \frac{{\partial p_{\text{f}} }}{\partial y}\frac{\partial y}{\partial t} = \frac{1}{cy}\frac{\partial y}{\partial t}$$
(42)

To express (∇p f)2, we should consider \(\nabla p_{\text{f}} = \partial p_{\text{f}} /\partial x\). Applying the chain rule for one dimension:

$$\frac{{\partial p_{\text{f}} }}{\partial x} = \frac{{\partial p_{\text{f}} }}{\partial y}\frac{\partial y}{\partial x}$$
(43)
$$\frac{{\partial p_{\text{f}} }}{\partial x} = \frac{{\partial p_{\text{f}} }}{\partial y}\nabla {\text{y }}$$
(44)

Substituting Eq. (40) in Eq. (44)

$$\left( {\nabla p_{\text{f}} } \right)^{2} = \frac{1}{{\left( {cy} \right)^{2} }}\left( {\nabla y} \right)^{2}$$
(45)

For the term ∇2 p f:

$$\nabla^{2} p_{\text{f}} = \left( {\frac{{\partial^{2} p_{\text{f}} }}{{\partial x^{2} }}} \right) = \frac{\partial }{\partial x}\left( {\frac{{\partial p_{\text{f}} }}{\partial x}} \right)$$
(46)

Substituting Eq. (43) into Eq. (46) gives:

$$\nabla^{2} p_{\text{f}} = \frac{\partial }{\partial x}\left[ {\left( {\frac{{\partial p_{f} }}{\partial y}} \right)\left( {\frac{\partial y}{\partial x}} \right)} \right]$$
(47)

Applying the derivative product and transposing terms give:

$$\nabla^{2} p_{\text{f}} = \frac{\partial }{\partial x}\left( {\frac{{\partial p_{\text{f}} }}{\partial y}} \right)\frac{\partial y}{\partial x} + \frac{\partial }{\partial x}\frac{\partial y}{\partial x}\left( {\frac{{\partial p_{\text{f}} }}{\partial y}} \right)$$
$$\nabla^{2} p_{\text{f}} = \frac{\partial }{\partial y}\left( {\frac{{\partial p_{\text{f}} }}{\partial x}} \right)\frac{\partial y}{\partial x} + \frac{{\partial^{2} y}}{{\partial x^{2} }}\left( {\frac{{\partial p_{\text{f}} }}{\partial y}} \right)$$
(48)

Substituting Eq. (43) into Eq. (48) gives:

$$\nabla^{2} p_{\text{f}} = \frac{\partial }{\partial y}\left[ {\left( {\frac{{\partial p_{\text{f}} }}{\partial y}} \right)\left( {\frac{\partial y}{\partial x}} \right)} \right]\left( {\frac{\partial y}{\partial x}} \right) + \frac{{\partial^{2} y}}{{\partial x^{2} }}\left( {\frac{{\partial p_{\text{f}} }}{\partial y}} \right)$$
$$\nabla^{2} p_{\text{f}} = \left( {\frac{{\partial^{2} p_{f} }}{{\partial y^{2} }}} \right)\left( {\frac{\partial y}{\partial x}} \right)^{2} + \frac{{\partial^{2} y}}{{\partial x^{2} }}\left( {\frac{{\partial p_{f} }}{\partial y}} \right)$$
(49)

Substituting Eqs. (40) and (41) into Eq. (49)

$$\nabla^{2} p_{\text{f}} = \frac{1}{cy}\left( {\nabla^{2} y} \right) - \frac{1}{{cy^{2} }}\left( {\nabla y} \right)^{2}$$
(50)

Substituting Eqs. (42), (45), and (50) into Eq. (20) gives:

$$\frac{1}{Dcy}\frac{\partial y}{\partial t} = \frac{1}{{c\left( y \right)^{2} }}\left( {\nabla y} \right)^{2} + \frac{1}{c\left( y \right)}\left( {\nabla^{2} y} \right) - \frac{1}{{c\left( y \right)^{2} }}\left( {\nabla y} \right)^{2}$$

Simplifying this equation, we obtained the linear diffusivity equation, Eq. (51).

$$\frac{1}{D}\frac{\partial y}{\partial t} = \nabla^{2} y$$
(51)

This equation is solved in Matthews and Russell (1967) for different conditions or cases: (1) constant rate, infinite reservoir; (2) constant rate, closed outer boundary; and (3) constant rate, constant pressure outer boundary case. Moreover, this type of equation is compared in papers (Chakrabarty et al. 1993; Odeh and Babu 1988) to validate these linear equations.

Equation (51) models fluid flow in a nonlinear diffusion process for a reservoir with extension fractures, without matrix–fracture transfer; this equation expressed for radial flow becomes:

$$\frac{{\partial^{2} y}}{{\partial r^{2} }} + \frac{1}{r}\frac{\partial y}{\partial r} = \frac{{\phi \mu^{*} c}}{k}\frac{\partial y}{\partial t}$$
(52)

Initial and boundary conditions for transformed Eq. (52) are:

  1. 1.

    y(r, 0) = y i at t = 0 for all r

  2. 2.

    \(\left( {r\partial y/\partial r} \right)_{{r_{\text{w}} }} = - 6q\mu /\pi ha^{2}\) for t > 0.

To develop the solution, this inner boundary condition is replaced by the line source condition:

$$\mathop {\lim }\limits_{r \to 0} 1/cy\left( {r\partial y/\partial r} \right)_{{r_{\text{w}} }} = - 6q\mu /\pi ha^{2} \,\,\,\,{\text{for }}\,\,t > 0$$
  1. 3.

    y(r, t) = y i as \(r \to \infty\) for all t.

The solution of Eq. (52) is:

$$y = \frac{{ y_{\rm i} }}{{1 - \frac{3q\mu c}{{\pi ha^{2} }}\left[ {E_{\rm i} \left( {\frac{{\phi \mu^{*} cr^{2} }}{4kt}} \right)} \right]}}$$
(53)

Substituting the initial transformation \(p_{\text{f}} = (1/c)\ln (cy)\) into Eq. (53):

$$p_{\text{f}} = \frac{1}{c}\ln \left( {\frac{{ y_{\text{i}} }}{{1 - \frac{3q\mu c}{{\pi ha^{2} }}\left[ {E_{\text{i}} \left( {\frac{{\phi \mu^{*} cr^{2} }}{4kt}} \right)} \right]}}} \right)$$
(54)

where E i represents the exponential integral function, for argument values <0.0025 [(4kt/ϕμ * cr 2) > 100]; the logarithmic approximation for the wellbore pressure is:

$$p_{\text{wf}} = \frac{1}{c}{ \ln }\left( {\frac{{ y_{\text{i}} }}{{1 - \frac{3q\mu c}{{\pi ha^{2} }}\left[ {\ln \left( {\frac{kt}{{\phi \mu^{*} cr_{\text{w}}^{2} }}} \right) + 0.80907} \right]}}} \right)$$
(55)

The resultant equations are given in “Appendix 2,” which also outlines the solution procedure.

A flow equation similar to Eq. (20) (nonlinear parabolic differential equation) was solved numerically using implicit finite difference for non-steady–steady flow, and the pressure distribution in fractured reservoirs (Yilmaz et al. 1994). In this paper, our solution is analytical, and it was developed applying Cole–Hopt and Boltzmann transformations.

Jump condition

A jump condition holds at a discontinuity or abrupt change. A composite media, £, is traditionally modeled using Darcy’s equation. At the interface between two media, the continuity of mass and momentum across the interface is required. So, a jump condition is needed.

Impermeable fractures have no jump of normal velocity while jumps of pressure are present in £; highly permeable fractures have jumps of normal velocity, with no pressure jump on £. In our study, the jump is in flow velocity. So, jump condition constrains the two states on either side of a discontinuity in conformance with conservation of momentum.

We adopt a Couette’s equation for the fractures, in which the volumetric flow rate q on £ satisfies Eqs. (4b) and (8a). Equation (10) satisfies the conservation of mass in the fracture.

Equation (56) guarantees the continuity of momentum:

$$\sigma^{\pounds} \left( {u,p} \right) = \sigma \left( u \right) - \alpha pI$$
(56)

It is necessary for the existence of a traction jump across a flat interface for dynamics of solid–solid phase transition. Here, I is the identity tensor, p is pressure, α is Biot constant associated with compressibility, σ(u) is the stress tensor, and u is velocity.

As already stated, fractures could be open or closed. If they are open fractures, a connected porous space will be observed in its width. However, a closed fracture does not have width. Then, \(n^{ - } = - n^{ + }\), which expresses velocity changes continuously from \(n^{ - } {\text{to }} - n^{ + }\) in porous media.

For any function g defined in £, the jump of g on £ in the direction of \(n^{ + }\):

$$\left[ g \right]_{{\text{\pounds}}}\,=\,g^{ + } - g^{ - }$$
(57)

The width w is the jump of \(u \bullet n^{ - }\) on £:

$$w = - \left[ u \right]_{\pounds}\,\bullet\,n^{ + } ,$$
(58)

where n is perpendicular to velocity (Chambat et al. 2014; Cermelli and Gurtin 1994).

Results and discussion

In this section, comparisons and a field example are presented to describe fluid flow in NFTRs.

Comparisons of linear wellbore pressure using Darcy and Couette flow: infinite reservoir. To examine the differences between Darcy and Couette flow without a square term, we used analytical models to calculate and compare wellbore pressures. Key data and reservoir parameters values employed are presented in Table 1.

Table 1 Data and parameters of a naturally fractured tectonic reservoir (N)
Full size table

Equation (25) is used, with its initial and boundary conditions, to describe the pressure distribution using Darcy law. Its solution is Eq. (29), namely Case 1A.

Equation (30) is used, with its initial and boundary conditions, to describe the pressure distribution using the Couette equation. Its solution is Eq. (32), namely Case 1B.

For the application of Darcy’s law in the diffusivity equation, and the Couette equation in the inner boundary condition, Case 1C, we employed Eq. (33) and its solution is given by Eq. (35).

For the above three cases, we used two types of geometry: slab of parallel fractures and single fracture, which are presented in Table 2.

Table 2 Application cases for NFTRs
Full size table

Figures 7 and 8 show pressure behavior for a single and slab fractures (smooth and rough) in a NFTR, considering Darcy and Couette equations for inner boundary condition and the diffusivity equation, for smooth fractures, wellbore fracture pressure decreased strongly because flow velocity is high, so that this assumption is not real for fractured reservoirs due to the friction effect between fracture walls.

Fig. 7
figure 7

Pressure behavior: Couette and Darcy flow for an equivalent fracture

Full size image
Fig. 8
figure 8

Pressure behavior: Couette and Darcy flow for slab fractures

Full size image

In this paper, we used Poiseuille’s law (Eq. 7), considering hydraulic conductivity, kh, and fracture aperture, a, for rough fractures. These parameters were determined using well testing for a NFCR (Singha and Al-Shammeli 2012). In contrast, we used Eq. (6) considering permeability, k, and fracture aperture, a, for smooth fractures.

Also, Figs. 7 and 8 display Case 1C, which were modeled using Couette equation as inner boundary condition, describing non-Darcy flow nearby the well; for a single fracture, the initial pressure drop is high and approaches a nearly constant pressure for long times. A similar behavior is also observed for a slab fracture.

When Darcy equation is used, pressure behavior in NFTR is quasi-constant. In other words, pressure drop is low, creating conservative results and overestimating fluid flow.

Conclusion

This paper has presented an analytical model for the description of fluid dynamics in naturally fractured tectonic reservoirs (essentially type I reservoirs, Nelson 2001). The models describe the pressure behavior and fluid flow in fractured non-stress-sensitive reservoirs, without matrix–fracture transfer.

From the preceding discussions and based on the material presented in this paper, the following conclusions can be made:

  1. 1.

    The analytic model presented an analysis of the single-phase flow equation for incompressible fluid, in non-stress sensitive naturally fractured tectonic reservoirs. This analysis showed the error in using the linear solution for NFTRs.

  2. 2.

    An analytical solution to quantify fluid dynamics in non-stress-sensitive NFTRs is proposed.

  3. 3.

    The nonlinear solution shows that for high flow rates there is a correction for the pressure and fluid flow, suggesting that the nonlinear term in Eq. (1) must be taken into account, to correctly describe oil flow through NFTRs due to non-Darcy laminar flow.

  4. 4.

    This study explains the phenomenon of high initial production rates, declining after a short period of time for fissured formations without interaction matrix–fractures.

  5. 5.

    Our results suggest that exact solution of Navier–Stokes equation, namely Couette’s flow, correctly predicts fluid flow in NFTRs.