© M. Abbasi et al., published by IFP Energies nouvelles, 2019

Nomenclature

A: Cross section area (m²)

a: Heterogenic function coefficient

c_f : Formation compressibility (psi⁻¹)

c_t : Total compressibility (psi⁻¹)

c_n : Series coefficient

F(x): Function

J_a : Bessel function of the first kind

k: Rock absolute permeability (md)

L: Matrix length (m)

P: Pressure (psi)

Q: Cumulative production (m³)

q: Wellbore flow rate (m³/day)

s: Laplace space

t: Time (s)

u,w: Change of the variable

x: Coordinate indicator

Y_a : Bessel function of the first kind

Greek

α: Permeability dependency coefficient (psi⁻¹)

μ: Fluid viscosity (cp)

∅: Rock porosity

ρ: Fluid density (kg L⁻³)

∂: Partial derivative symbol

σ: Shape factor (L⁻²)

Subscripts

D: Dimensionless

i: Initial

f: Fracture

m: Matrix

w: Well

1 Introduction

Fractured reservoirs contain large percentage of the worldwide oil. There are two different matrix and fracture systems with different physical characteristics inside such oil reserves. The mechanisms of production of this type of reservoirs are very different from conventional reservoirs (Rangel-German et al., 2010). Considering the juxtaposition of two matrix and fracture systems together (with considerable difference in size), the simulation of fractured reservoirs using single porosity models is very time-consuming and even the simulation may not converge. Hence, modeling of this type of reservoirs via dual porosity models with the definition of the fluid transfer coefficient between the matrix block and the fracture is much more common (Abbasi et al., 2018b).

The dual porosity model was introduced for the first time by Barenblatt et al. (1960). In this model, the fractured reservoir consists of two systems: the matrix block that stores oil and the fracture that is responsible for conducting the flow of fluid to the well. Then, Warren and Root (1963) interpreted the data of this type of reservoirs by developing the relations of the Barenblatt et al. (1960) specifically in petroleum engineering.

The introduction of dual-medium dates back to 1960s (Barenblatt et al., 1960; Warren and Root, 1963) to facilitate modeling the fluid flow in vugular or naturally fractured reservoirs which consist of two rather different media including fracture (f) and matrix (m) systems. The governing flow equations in such a system may be represented by the following two major equations assuming slightly compressible fluid:

(1)

(2)which describe the fluid flow in fracture and matrix. Barenblatt and Zheltov presumed that fracture porosity does not contribute to the system. In the dual porosity model introduced by Warren and Root, the only contributor to the fracture system is the matrix system. Thus, the second equation is reformed as:

(3)

Equation (3) presents the matrix-fracture volumetric flow exchange per matrix bulk volume, and contains an exchange coefficient called shape factor (σ). The previous model presumes a PSS matrix-fracture shape factor. σ is a fixed coefficient and is proportional to inverse of square of matrix dimensions. Warren and Root formulated a mathematical expression to express the shape factor for the fractured reservoirs consisting of cubical matrix blocks separated by orthogonal fractures uniformly.

Kazemi (1969) used numerical simulator to define the shape factor between matrix and fracture system. Gilman and Kazemi (1983) developed a simulator in heterogeneous fractured reservoirs using their shape factor. Zimmerman et al. (1993) used the Fourier series to present an analytic relation for the shape factor of different geometric shapes. Lim and Aziz (1995) derived matrix–fracture shape factors for single-phase flow by solving the diffusivity equation in a cylindrical and spherical geometry, respectively. They also expressed that the matrix geometry has a great impact on the matrix-fracture transfer shape factor. Bourbiaux et al. (1999) considered the fluid transfer coefficient between the matrix and the fracture in two dimensions and pseudo steady state conditions by presenting the expression of the shape factor (or exchange coefficient) in terms of the dimensions of the block of matrix. Noetinger and Estebenet (2000) calculated matrix-fracture transient and steady-state exchange terms in a naturally fractured porous medium based on Continuous Time Random Walk (CTRW) methods. Noetinger et al. (2001a) presented a novel formation of CTRW algorithm can depict a direct determination of transient exchange function vastly used in the interpretation of well testing analysis based on probability density function. In Noetinger et al. (2001b) study, the identity based upon two definitions of steady state exchange coefficient due to moderately disparate physical methodologies. One definition is to consider a steady state closure system, while the other being intrinsically transient. The objective is to prove that both methods essentially give the same well-defined value of coefficient α. In Landereau et al. (2001) paper, the behavior of large-scale characteristics of 2D highly fractured systems is put forth via deterministic up-scaling technique. Firstly, large-scale fracture permeability tensor is examined, which is affiliated with the interconnection of fracture systems. Secondly, the dependency of exchange on matrix block geometry is investigated. At last, the proposed methodology is compared with several methods published in the literature.

Germán (2002) presented the experimental and analytical data to study the concept of the fluid transfer coefficient between the matrix and the fracture. Rangel-German and Kovscek (2006) examined the mass transfer between the matrix and the fracture system attributed to the fluid flow due to the difference in pressure and imbibition due to capillary force. Hassanzadeh and Pooladi-Darvish (2006) studied the effect of the boundary conditions of the matrix block on the shape factor between the matrix block and the fracture in different geometries of the flow. Hassanzadeh et al. (2009) studied the flow modeling in a radial fractured reservoir with constant rate condition and calculated the transient shape factor in the dual porosity model by combining the flow-in-fracture and matrix flow equations. Bourbiaux and Ding (2016) provided a technique to assess the production from unconventional low-permeability reservoirs considering the transient effects as well as transfer functions embedded in double porosity simulators. Abbasi et al. (2018a) investigated the effect of block size distribution on the fluid flow and shape factor between the matrix and fracture by presenting an analytical model. In most of the abovementioned works, simplifying assumptions were made for developing dimensionless shape factor. Ding et al. (2018) proposed an Embedded Discrete Fracture Model – Multiple INteracting Continua (EDFM-MINC) methodology to model shale reservoirs by specifying equivalent porosity and permeability for the fractures with low-conductivity property. In addition, a number of methods are introduced to improve transmissibility calculations among several media of an EDFM by taking into account the pressure distribution using an integral method. In Farah et al. (2019) paper, a novel doable hybrid method based upon a MINC proximity function is put forth to predict flow behavior in fractured reservoirs of low permeability type. The proposed method has the advantages of both discrete fracture model and multi-continuum for the purpose of simulating unconventional fractured reservoirs. It is mainly suitable for interference test analysis, large reservoir modeling, and/or when lack of data exacerbate the possibility of sensitivity studies.

In line with previous studies, fluid flow modeling in porous media based on dual porosity concept involve assumptions that distance them from real prevailing conditions in reservoirs. Such assumptions include ignoring quadratic pressure gradient term, homogeneous porous media, and constant permeability during the reservoir life. In reservoirs with low permeability, and high viscous compressible fluids, the quadratic pressure gradient term affects the results and produces erroneous predictions in well testing particularly within the near wellbore regions where the pressure drop is high. Moreover, the assumption of matrix-fracture permeability makes the dual porosity based simulators predict smaller pressure drop, since in real reservoir conditions, matrix and especially fracture permeability reduces as oil production occurs. This paper investigates the effect of reducing matrix permeability on the pressure drop.

In this paper, considering more realistic assumptions about the nature of the fractured reservoirs such as the heterogeneity of the matrix block, the quadratic gradient pressure in the equation of diffusion and the change of rock properties by changing the pressure, the analytical modeling of the fluid flow from the matrix to the fracture is discussed. The effect of the parameters expressed on the shape factor is examined and compared with the previous models. In order to validate the results of the proposed analytical model, results of the new method are compared with the results of the commercial simulator.

2 Basic model

In this section, the diffusivity equations in the matrix block are discussed. Additionally, considering the fracture as a boundary condition and taking into account the assumptions close to the nature of fractured reservoirs, the governing equations for fluid flow are introduced (Fig. 1). The diffusivity equations within the matrix block are expressed in terms of three parameters of quadratic pressure gradient, heterogeneous porous medium and pressure-dependent rock properties.

2.1 Quadratic pressure gradient term

One of the simplification assumptions that has been considered so far is the underestimation of quadratic pressure gradient term in which case the pressure in the low porosity porous medium modeling such as the matrix block causes an error and will affect the shape factor. The diffusivity equation in the matrix block has been derived in Appendix which is similar to the derivation process reported by Stewart (Abbasi et al., 2017; Ahmed and McKinney, 2011; Stewart, 2011; Zimmerman, 2002).

(4)

To make this mathematical model more tractable and easier to understand, the following assumptions are given:

• homogenous and isotropic porous media,

• constant fluid properties,

• flow of slightly compressible fluid,

• isothermal conditions.

2.2 Pressure-dependent matrix properties

By decreasing the reservoir pressure with respect to the compressibility of the matrix block, the static properties of the matrix block such as porosity and permeability are reduced. Applying the pressure-dependent matrix properties, the diffusivity equation in the matrix block is expressed as:

(5)

2.3 Heterogeneous porous medium

The flow equations in fractured reservoirs are based on the homogeneity of the porous medium, while in view of the nature of oil reservoirs and sedimentation of formations in different geological periods, the porous medium is heterogeneous and the permeability of the matrix block varies along the block. The flow equation in the heterogeneous matrix block is expressed as:

(6)

3 Analytical model

In this section, considering the fracture as the boundary condition, the analytical solution of the diffusivity equation is presented taking into account the quadratic pressure gradient, the heterogeneous permeability and the pressure-dependent properties of the rock.

3.1 Quadratic pressure gradient term

Equation (4) is a partial differential equation and related initial and boundary conditions are defined as follows:

(7a)

(7b)

(7c)

In order to simplify the expressed equations, the following dimensionless parameters are used:

(8a)

(8b)

(8c)

Using dimensionless parameters, equations (4) and (7a)–(7c) are expressed as:

(9a)

(9b)

(9c)

(9d)where η is defined as:

Using the change of the following variable:

(10)Equations (9a)–(9d) are expressed as:

(11a)

(11b)

(11c)

(11d)

Using Laplace transform, we have:

(12a)

(12b)

(12c)

The general solution of equations (12a)–(12c) is as follows:

(13)

Using the boundary conditions (12b) and (12c), equation (13) will be:

(14)

Using equation (10) in Laplace space, we have:

(15)

3.2 Pressure-dependent matrix properties

In order to solve the partial derivative partial differential equation of fluid in a porous medium with pressure-dependent properties (Eq. (5)) in terms of the compressibility of the matrix block, porosity- and pressure-dependent permeability can be expressed as follows (Wang et al., 2018):

(16a)

(16b)

Given the porosity and permeability functions, equation (5) is expressed as:

(17)

Using the following change of variable:

(18)

Equation (17) is expressed as follows:

(19)

Given the initial and boundary conditions equations (7a)–(7c) and the change of variable equation (18), the initial and boundary conditions of equation (19) are as follows:

(20a)

(20b)

(20c)Using Laplace transform:

(21a)

(21b)

(21c)

The general solution of equation (21a) is as follows:

(22)

Using the boundary conditions (21b) and (21c), equation (22) will be:

(23)

Using equation (18) in Laplace space, we have:

(24)

3.3 Heterogeneous porous medium

In order to solve the partial differential equation of the fluid flow in a heterogeneous porous medium (Eq. (6)), two linear and nonlinear functions of permeability changes relative to the location were considered. The permeability function in equation (3) can be expressed in two ways:

(25a)

(25b)

3.3.1 Nonlinear function

Using the dimensionless definitions of equations (8a)–(8c) and the nonlinear function of permeability changes based on the location, equation (6) is as follows:

(26)

The initial and boundary conditions of equation (26) will be similar to equations (9b)–(9d). Given the following change of variable:

(27)

Equation (26) and its initial and boundary conditions are expressed as follows:

(28a)

(28b)

(28c)

(28d)

Then, using the following change of variable:

(29)

Equations (28a)–(28d) are expressed as:

(30a)

(30b)

(30c)

(30d)

Using Laplace transform we have:

(31a)

(31b)

(31c)

The general answer to equation (31a) is as follows:

(32)

Using the boundary conditions (31b) and (31c), equation (32) will be:

(33)and will further be expressed in terms of dimensionless location as:

(34)

3.3.2 Linear function

If the permeability variation of the porous medium is linear, using the dimensionless definitions of equations (8a)–(8c), the partial derivative of equation (6) is expressed as follows (Kuchuk et al., 1996):

(35)

The initial and boundary conditions of equation (35) will be similar to equations (9b)–(9d). By changing the following variable:

(36)

Equation (35) and its initial and boundary conditions are expressed as follows:

(37)

Equation (37) is similar to the radial flow equation of fluid in a porous medium. Given the applied change of variable, the initial and boundary conditions of equation (37) is equal to:

(38a)

(38b)

(38c)

Using the Laplace transform of equations (38a)–(38c) and its boundary conditions are expressed as follows:

(39a)

(39b)

(39c)

The general answer to equation (39a) will be as follows:

(40)

Using the boundary conditions (39b) and (39c), equation (40) is expressed as follows:

(41)

Equation (41) can be expressed in terms of x_D using equation (36).

4 Result and discussion

In this section, using the analytical solutions derived in the previous section, the effect of the parameters of the quadratic pressure gradient, the heterogeneous porous medium and the pressure-dependent rock properties on the oil rate and the shape factor between the matrix and the fracture are examined. To determine the amount of fluid transferred between the matrix block and the fracture, the dimensionless shape factor with the following mathematical description is utilized (Hassanzadeh and Pooladi-Darvish, 2006):

(42)

4.1 The effect of quadratic pressure gradient term

Figure 2 shows the dimensionless flow rate plot from the matrix block to the fracture. As seen in the figure, by increasing the quadratic pressure gradient term η (in other words, considering the effect of the quadratic pressure gradient), the pressure drop in the matrix block decreases and the fluid transfer rate increases in a constant time (i.e., higher rate for larger value of η).

Figure 3 shows the dimensionless shape factor between the matrix block and the fracture. Considering the quadratic pressure gradient in the diffusion equation, the rate of fluid transfer between the matrix block and the fracture and the transient dimensionless shape factor increases.

Considering the expression of the quadratic pressure gradient, the pressure drop in the matrix block also decreases, which delays the pressure drop to reach the end of the matrix block and increases the transient time period.

It can also be said that by increasing the pressure dependency of the matrix rock and consequently increasing the flux of the matrix block, the effect of this expression increases.

In Figure 4, the dimensionless pressure of no-flow limit within the matrix block (x = L/2) and dimensionless shape factor is shown. As shown in the figure, with the arrival of the pressure drop to the end of the matrix block, three time regions can be defined: (A) the time interval when the pressure drop has not yet reached the end of the matrix block and the dimensionless shape factor changes linearly versus time; (B) the pressure drop reaches the end of the matrix block, and pressure of no-flow limit within the matrix block is decreasing and the dimensionless shape factor curvature is reduced; (C) pressure of no-flow limit within the matrix block is equal to the fracture pressure and no fluid transitions occur between the matrix and the fracture. In this case, the dimensionless shape factor reaches its lowest value and is set to 9.87. Lim and Aziz (1995) have derived the same stabilized value. With regard to the pressure gradient quadratic term, the time of arrival of the pressure drop is delayed, which causes the time period of the transient region to increase.

4.2 The effect of pressure-dependent matrix properties

Figure 5 shows the permeability of the matrix block at t = 100 s. As seen in the figure, in the areas close to the fracture (x = 0), the pressure drop is higher, the permeability decrease is relatively higher, and the permeability of the matrix block increases in the regions far from the fracture. Figure 6 shows the oil flow rate versus time. As shown in the figure, taking into account the permeability dependence of the matrix block on the pressure, the permeability of the matrix block decreases with the production of oil, which causes the amount of oil input from the matrix block to the fracture to be less than the constant permeability model at a fixed time.

Figure 7 shows the dimensionless shape factor between the matrix block and the fracture. In the model of fluid flow in a porous medium, assuming a constant permeability, the amount of oil outflow from the matrix block is greater than other models, which corresponds to a higher value of the dimensionless shape factor during the transient period of matrix-fracture transfer. Taking into account the pressure-dependent permeability of the matrix, the pressure drop reaches the end of the matrix block with delay, and as shown in Figure 7, the dimensionless shape factor increases at the transient period, and after the arrival of the pressure drop to the end, the dimensionless shape factor between the matrix and the fracture is fixed.

4.3 The effect of heterogeneity in porous medium

Figure 8 shows the permeability of the matrix block versus location. The parameter a in linear and quadratic functions, respectively, indicates the line slope and curvature of the function, so that with increasing a, the variation of permeability relative to the location increases.

The defined shape factor in equation (42) is based on dimensionless shape factor explained in equation (8c) considering a porous medium with fixed permeability. In order to compare shape factors associated with different permeability, one can make use of the corrected form of equation (42) with average permeability.

(43)

Figure 9 shows the dimensionless shape factor between the matrix and the fracture for the three different permeation models. The graph of dimensionless shape factor versus dimensionless time has been plotted considering permeability in homogenous porous media. As depicted in the figure, as matrix permeability reduces, the dimensionless shape factor increases within the transient region. If the corrected graph of dimensionless shape factor versus dimensionless time is plotted based on the mean permeability of the porous medium, all three graphs will match.

5 Validation

5.1 Analytical model

A commercial numerical simulator was used to validate the analytical model. The schematic model of the single porosity model is shown in Figure 10. The fine single porosity model has dimensions of 948 ft × 948 ft × 32/8 ft and the number of grids is 31 × 31 × 1. In order to verify the one-dimensional flow, the matrix permeability was considered zero in the y and z directions, and the fracture system was considered as a very large grid having a permeability of 10 000 md.

The parameters used to simulate the single porosity model are expressed in Table 1. In order to investigate the effect of permeability heterogeneity, three models different in permeability were considered. As shown in Figure 11, the permeability of the matrix block is plotted in three models: homogeneous, linear and the quadratic permeability function is plotted according to the location (grid permeability decreases as we go away from the fracture location).

Figure 12 shows the matrix block average pressure of the analytical model and numerical simulator versus time. As it is shown in the figure, the matrix block average pressure by decreasing the permeability of the matrix block with two linear and quadratic heterogeneity functions is greater than the homogeneous model. Based on the two heterogeneous models, as we step away from the fracture, the permeability decreases; at early times when the pressure drop occurs in the grids near the fracture, the results of all three models are the same, and with increasing time and the arrival of the pressure drop on the upper grids, this difference increases. It can also be said that in heterogeneous models, especially the quadratic heterogeneity model, the permeability at the upper part of the model is negligible, which results in small depletion and consequently, the pressure remains almost constant and equal to the initial pressure. The average block pressure of the matrix will be higher than the homogeneous model. The results of the analytical model show very good agreement with the results of the numerical simulation.

5.2 Time-dependent shape factor

To validate the dimensionless shape factor between the matrix and the fracture in the dual porosity model of the numerical software, the decreasing trend in dimensionless shape factor obtained in the homogeneous model, which consists of two transient and pseudo steady state parts, was used and the results were compared with the single porosity model. The numerical software uses Kazemi model (Gilman and Kazemi, 1983) for the shape factor, which considers a constant value, and it can be used in each time step taking into account a coefficient for the shape factor to apply the transient shape factor. In order to calculate the factor associated with dimensionless shape factor in first time step, the average dimensionless shape factor is assessed from the graphs provided in this paper in first time step, and is further divided by the dimensionless shape factor in the pseudo state flow regime. In addition, the average dimensionless shape factor in the second time step () is divided by the already calculated first time step dimensionless shape factor () to assess the second time step dimensionless shaper. As put forth in Table 2, other factors related to dimensionless shape factor are calculable in each time step. It is worth noting that as time passes and appearance of pseudo steady state region, the graph of dimensionless shape factor becomes a straight line, which implies that the value due to division of average dimensionless shape factor in each time step by the previous time step becomes more unity.

As shown in Figure 13, considering the coefficients defined for the dimensionless shape factor at each time step in Table 3, the recovery matrix oil recovery of the dual porosity model (with transient shape factor) matches well the single porosity model (model part a).

6 Conclusion

In the dual porosity reservoirs, shape factor is used to express the fluid transfer rate from the matrix to the fracture system. The previous models for shape factors focused more on the pseudo steady state part of the flow and expressed a constant value for the shape factor; however, the fact that the flow rate of the matrix block to the fracture in the pseudo steady state part is very small and most of fluid transfer between the matrix block and the fracture occurs in the transient region has not been investigated.

In the present study, considering the quadratic pressure gradient, the heterogeneous matrix and the pressure-dependent rock properties, an analytical solution of the flow of the fluid from the matrix to the fracture was derived. Considering the quadratic pressure gradient, the pressure drop in the matrix block decreases and this causes the shape factor between the matrix and the fracture in the transient region to increase as compared to the previous models. The effect of this parameter increases with increasing the compressibility of the rock and fluid. Considering the permeability as a function of location and pressure, the amount of fluid transfer between the matrix and the fractures changes, and consequently the dimensionless shape factor value changes. Based on the analytical relationships expressed for the dimensionless shape factor, new calculated values can be used in commercial simulators to improve the dual porosity model.

Diffusivity equation considering quadratic pressure gradient

The general form conservation of mass equation in Cartesian coordinate is (Abbasi et al., 2017):

(A1)

Diffusion equation which is certain form of partial differential equation governs the transient flow through porous media. To achieve this equation, Darcy’s law, mass conservation equation, and an equation which explicates the fluid inside rock is combined together. Considering the right hand side of equation (A1), and utilizing the chain rule of differentiation:Darcy’s law for fluid flow in porous media is as follows:

(A2)

Now look at the left-hand side of equation (A1). The flux q is given by Darcy’s law:

(A3)

By combination of above equations, general form of diffusivity equation considering the quadratic pressure gradient is obtained as follows:

(A4)

References

All Tables

All Figures