Quasi-Analytical Model of the Transient Behavior Pressure in an Oil Reservoir Made Up of Three Porous Media Considering the Fractional Time Derivative

Abstract

The present work proposes a new model to capture high heterogeneity of single phase flow in naturally fractured vuggy reservoirs. The model considers a three porous media reservoir; namely, fractured system, vugular system and matrix; the case of an infinite reservoir is considered in a full-penetrating wellbore. Furthermore, the model relaxes classic hypotheses considering that matrix permeability has a significant impact on the pressure deficit from the wellbore, reaching the triple permeability and triple porosity model wich allows the wellbore to be fed by all the porous media and not exclusively by the fractured system; where it is considered a pseudostable interporous flow. In addition, it is considered the anomalous flow phenomenon from the pressure of each independent porous medium and as a whole, through the temporal fractional derivative of Caputo type; the resulting phenomenon is studied for orders in the fractional derivatives in (0, 2), known as superdiffusive and subdiffusive phenomena. Synthetic results highlight the effect of anomalous flows throughout the entire transient behavior considering a significant permeability in the matrix and it is contrasted with the effect of an almost negligible matrix permeability. The model is solved analytically in the Laplace space, incorporating the Tartaglia–Cardano equations.

1. Introduction

In order to create an efficient plan to exploit a reservoir, proper design facilities at the bottom of a well and in the surface and get a more realistic overview of the conditions of wellbore and of the reservoir that help one to make better decisions to preserve the injectivity and later secondary recovery of the fluid, it is necessary to correct reservoir modeling that allows the cover of the heterogeneity of the fluid. Therefore, pressure tests are one of the main tools to determine reservoir characteristics and at the same time determine the main factors that influence hydrocarbon production.

As is well known, Warren and Root [1] were among the pioneers in work with naturally fractured reservoirs and started modeling non-homogeneous reservoirs, proposed a double porosity model and thought of the reservoir as a set of rectangular parallelopipeds (which represents the matrix rock), while the secondary porosity, the fractures, form a continuous orthogonal system uniformly spaced by the parallelepipeds; they considered the fluid in the matrix as being in a quasi-stable state flow while the fluid in the fractures is unstable.

Adams et al. [2] analyzed the pressure tests in the buildup phase and emphasized two main slopes in it, the first one related to the permeability of the matrix and the second one related to the permeability of the fracture-matrix system.

Kazemi [3], on the other hand, generalized Warren and Root’s work by assuming that the fluid throughout the reservoir is unstable; in addition, Kazemi considered in his model that the fractures are all horizontal and uniformly distributed.

Moench [4] studied the works of Warren and Kazemi and developed a theory that incorporates the term of hydraulic skin present in the matrix-fracture interface, which allows one to generalize the assumptions of pseudo-stable Warren–Root [1] fluid and fully transient flow of Kazemi [3]. Moreover, the skin is defined to characterize the condition of the well and the degree of connectivity between the well and the reservoir, the result being that, for a positive skin $s>0$, damage is considered in the well where the greater the damage, the greater the value of the skin; while a negative skin $s<0$ represents a stimulated well [5]. Another important concept is the well storage effect, which occurs when fluid is compressed into a completely filled well; characterized by causing a straight line behavior with slope 1 in the times immediately following the start of production [6].

Several observations in pressure well-tests showed that double porosity models were insufficient to explain and reproduce these behaviors; so Abdassah and Ershaghi [7] developed a triple porosity model formed by two types of matrix separated by horizontal layers of fractures and showed that triple porosity systems are characterized by various anomalous slope changes generated by the respective contribution of each matrix block.

Due to Bourdet [8], who developed a new method based on the analysis of the time rate of pressure change, diagnosing pressure-transient behavior was simplified to identify characteristic regimes that produce straight lines as Ehlig emphasized [9].

Analyzing the behaviors in naturally fractured vuggy carbonate reservoirs, Camacho-Velázquez et al. [10] and Fuentes-Cruz et al. [11] proposed a triple porosity and double permeability model for fully and partially penetrating wells. They considered the porous medium as a continuum constituted by fractures, matrix and vugs (where it does not distinguish between channels, caverns and vugs); they distinguished that the presence of vugs affects the behavior of pressure curves. In their model, they allowed an interaction between each of the porous media using the pseudo-stable flow approximation.

Following this approach, where the porous medium is constituted by fractures, matrix and vugs, Wu et al. [12] proposed a model where the flow is radial and only the fractures feed the producing well.

Thus, different attempts to capture the heterogeneity in fractured reservoirs were made, varying the considerations of flow, geometry or idealization of types of porous media that make up the reservoir. Al-Ahmadi et al. [13] proposed a triple porosity model considering the matrix rock and two different types of fractures (little permeable microfractures and more permeable macro fractures) for pseudo-stable and transient flows periods for linear and radial geometries.

Nie et al. [14] addressed the triple porosity model considering the shape of vugs as spherical, and therefore the interporous flow between fracture and vugs is unsteady, but the interporosity flow between matrix and vugs was ignored.

Some variants made in order to make the existing models more robust were made, such as Camacho-Velazquez et al. [15], who resumed their proposed model for fully and partially penetrated deposits incorporating information from pressure tests, imagelogs, cores, well logs, among others. In this work, they concluded that vertical communication between vugs may be more important than their horizontal communication in partially penetrated wells. Or through optimization methods, like Gómez et al. [16], who implemented the tunneling method to find the parameters in the triple porosity double permeability model and given the problem of multiple minimums, they suggested the use of additional information to specify the physically correct solution.

Considering an anomalous diffusion of the fluid throughout the whole reservoir, Metzler et al. [17] showed that this behavior can be modeled through the inclusion of fractional derivatives, allowing the introduction of a memory effect.

In spite of the work developed by Metzler, until that moment, understanding of the geometric and physical interpretation of the fractional operators was minimal; thus, Moshrefi-Torbati and Hammond [18] concluded that fractional operators can be grouped as filters with partial memory varying between the extremes of complete memory and non-memory, depending on the order of the fractional derivative.

Caputo [19], by studying the variation in size in the pores through which the fluid moves, modified Darcy’s law by introducing the memory effect through the inclusion of the fractional time derivative; however, in order to maintain a dimensional balance, he modified the definition of permeability by introducing his term pseudopermeability.

Considering the modified Darcy’s law and combining it with the continuity equation, Raghavan [20,21] developed a model to understand fluid transport in tight rocks.

From another point of view, theoretical advances have been made, generalizing areas such as the dimension of the environment in which the reservoir is embedded, that is, considering that the space through which the fluid moves can approximate as a fractal. Razminia et al. [22] combined fractal geometry and fractional calculus to analyze the pressure in double porosity reservoirs, by integrating the fractional derivative in a fractal system considering the order of the fractional derivative congruent for both porous systems; more recently [23], they analyzed pressure in partially penetrated wells using Green’s formula; in both cases, the dimensional balance is ignored.

Posadas-Mondragón and Camacho-Velázquez [24] used this approach to solve the model for a partially penetrated well, considering that the distribution and connectivity of the fracture network as a fractal and the reservoir as a single porous medium.

On the other hand, Tian et al. [25] investigated the transient pressure behavior in a horizontal well in coalbed methane reservoirs with multiple fracture wings by identifying the characteristic flow periods. Meanwhile, Gao et al. [26] addressed the problem of fast pressure drop and gas well abandonment caused by well interference; divided the region around the MFHW into three regions and solved the model using finite volumes and concluded that the IWCM can effectively control the resources within the reservoir around the well.

Chang et al. [27], starting from the fact that integer derivatives cannot reliably quantify fluid flow with time memory and/or spatial path dependency, proposed and evaluated a Darcy’s law model with a spatial fractional derivative in the sense of Riemann-Liouville; showed how their model can felicitously depict flow properties; explored subdiffusive flows and saw how these flows could represent porous media with small pore sizes or dead pores.

Gu et al. [28] developed a double porosity model for multiple horizontal wells in a tight gas reservoir based on the Warren and Root model, incorporating fractal fracture network and anomalous behavior of diffusion processes through fractional calculus and explored the subdiffusive phenomenon.

Jiang et al. [29] considered an anomalous diffusion in multi-fractured horizontal wells in coalbed methane reservoirs, building a model that integrates fractal theory and fractional calculus and considered the stress sensitivity of the cleat systems; identified the characteristic flow regimes and summarized the influence laws of the parameters.

Thus, in this work, a further effort is made to capture the high heterogeneity in naturally fractured vugular reservoirs using fractional calculus as a tool. The idealization is resumed that the deposit is constituted by three independent and interconnected porous media in which the memory effect phenomenon is studied; however, it is allowed that the anomalous behavior can be proper and independent of each porous medium where tight matrix restriction is eliminated and therefore the possibility of a significant permeability in the matrix is allowed.

The work consists of Section 2, where the deduction of the dimensional fractional model is shown; Section 3, where the dimensionless model, called the triple permeability and triple porosity fractional model, is presented and solved, the model preserves a dimensional balance and the superdiffusive and subdiffusive phenomena are distinguished for each porous medium; Section 4, where synthetic results are shown, showing the impact of the anomalous flow in each phase of the transitory phenomenon and it is contrasted with the effect generated by a significant permeability in the matrix, called triple permeability, together with the impact generated by the inclusion of the constant that allows maintaining the dimensional balance; Section 5, which summarizes the conclusions reached in this work. Finally, Appendix A shows the asymptotic behavior in short times in the pressure deficit.

2. Deduction of Fractional Triple Permeability and Triple Porosity Model

Assume that the reservoir is made up of three porous media: fractured media, vugular media and matrix. Let $\varphi$ the porosity of the whole reservoir, ${\varphi }_f$ the porosity of the fractured system, ${\varphi }_v$ the porosity of the vuggy system and ${\varphi }_m$ the porosity of the matrix; then

$$\varphi ={\varphi }_m+{\varphi }_f+{\varphi }_v.$$

(1)

Let $\mathbf{q}$ be the Darcy flow of the whole reservoir, ${\mathbf{q}}_{\mathbf{f}}$ the Darcy flow of the fractured system, ${\mathbf{q}}_{\mathbf{v}}$ the Darcy flow of the vuggy system and ${\mathbf{q}}_{\mathbf{m}}$ the Darcy flow of the matrix, then

$$\mathbf{q}={\mathbf{q}}_{\mathbf{m}}+{\mathbf{q}}_{\mathbf{f}}+{\mathbf{q}}_{\mathbf{v}}.$$

(2)

Now, the generic continuity equation is

$$\frac{\partial \left(\rho \theta \right)}{\partial t}+\nabla ·\left(\rho \mathbf{q}\right)=\rho {\rm Y}$$

(3)

where $\rho$ is density, ${\rm Y}$ is a term of source or sink and $\theta$ is the volumetric fluid content with $0<\theta <\varphi$; note that, when the porous medium is saturated with fluid $\theta =\varphi$. The Darcy’s law is

$$\mathbf{q}=-\frac{1}{\mu }k\left(p\right)\left(\nabla p-\rho g\nabla D\right);$$

(4)

where g is the gravitational acceleration and D is depth as a function of spatial coordinates.

Then, the general flow transfer equation results from the combination of Equations (3) and (4).

Therefore, assuming that the reservoir is saturated with fluid, the general flow transfer equation for each porous medium is

$$\begin{array}{c}\hfill \frac{\partial \left(\rho {\varphi }_m\right)}{\partial t}+\nabla ·\left(\rho {\mathbf{q}}_{\mathbf{m}}\right)=\rho {{\rm Y}}_m,{\mathbf{q}}_{\mathbf{m}}=-\frac{1}{\mu }{k}_m\left(p_m\right)\left(\nabla p_m-\rho g\nabla D\right);\end{array}$$

(5)

$$\begin{array}{c}\hfill \frac{\partial \left(\rho {\varphi }_f\right)}{\partial t}+\nabla ·\left(\rho {\mathbf{q}}_{\mathbf{f}}\right)=\rho {{\rm Y}}_f,{\mathbf{q}}_{\mathbf{f}}=-\frac{1}{\mu }{k}_f\left(p_f\right)\left(\nabla p_f-\rho g\nabla D\right);\end{array}$$

(6)

$$\begin{array}{c}\hfill \frac{\partial \left(\rho {\varphi }_v\right)}{\partial t}+\nabla ·\left(\rho {\mathbf{q}}_{\mathbf{v}}\right)=\rho {{\rm Y}}_v,{\mathbf{q}}_{\mathbf{v}}=-\frac{1}{\mu }{k}_v\left(p_v\right)\left(\nabla p_v-\rho g\nabla D\right).\end{array}$$

(7)

To get back the flow transfer equation for the whole reservoir from Equations (5)–(7), it must be satisfied that

$${\rm Y}={{\rm Y}}_m+{{\rm Y}}_f+{{\rm Y}}_v.$$

(8)

Assuming that in the reservoir there is no gain or loss of fluid, then ${\rm Y}=0$ and consequently

$$\begin{array}{c}\hfill {{\rm Y}}_m={{\rm Y}}_{mf}+{{\rm Y}}_{mv};\end{array}$$

(9)

$$\begin{array}{c}\hfill {{\rm Y}}_f=-{{\rm Y}}_{mf}+{{\rm Y}}_{vf};\end{array}$$

(10)

$$\begin{array}{c}\hfill {{\rm Y}}_v=-{{\rm Y}}_{mv}-{{\rm Y}}_{vf}.\end{array}$$

(11)

Finally, assuming pseudostable interporous flow

$$\begin{array}{c}\hfill {{\rm Y}}_{mf}=a_{mf}\left(p_f-p_m\right);\end{array}$$

(12)

$$\begin{array}{c}\hfill {{\rm Y}}_{mv}=a_{mv}\left(p_v-p_m\right);\end{array}$$

(13)

$$\begin{array}{c}\hfill {{\rm Y}}_{vf}=a_{vf}\left(p_v-p_f\right).\end{array}$$

(14)

Hence, with $D=z$, the triple permeability and triple porosity model is

$$\begin{array}{c}\hfill {\varphi }_m{c}_m\frac{\partial p_m}{\partial t}=\nabla \left(\frac{k_m}{\mu }\nabla p_m\right)+a_{mf}\left(p_f-p_m\right)+a_{mv}\left(p_v-p_m\right),c_m=\frac{1}{{\varphi }_m}\frac{\partial {\varphi }_m}{\partial p_m};\end{array}$$

(15)

$$\begin{array}{c}\hfill {\varphi }_f{c}_f\frac{\partial p_f}{\partial t}=\nabla \left(\frac{k_f}{\mu }\nabla p_f\right)-a_{mf}\left(p_f-p_m\right)+a_{vf}\left(p_v-p_f\right),c_f=\frac{1}{{\varphi }_f}\frac{\partial {\varphi }_f}{\partial p_f};\end{array}$$

(16)

$$\begin{array}{c}\hfill {\varphi }_v{c}_v\frac{\partial p_v}{\partial t}=\nabla \left(\frac{k_v}{\mu }\nabla p_v\right)-a_{mv}\left(p_v-p_m\right)-a_{vf}\left(p_v-p_f\right),c_v=\frac{1}{{\varphi }_v}\frac{\partial {\varphi }_v}{\partial p_v}.\end{array}$$

(17)

Including the memory phenomenon in the dynamics of fluid movement, it is considered that the flow in the reservoir may be anomalous, so that the fractional derivative is included in the time derivative in the flow equation of each porous medium, obtaining the following system of partial differential equations.

$$\begin{array}{c}\hfill {\varphi }_m{c}_m{\tau }^{{\alpha }_m-1}\frac{{\partial }^{{\alpha }_m}{p}_m}{\partial t^{{\alpha }_m}}=\nabla \left(\frac{k_m}{\mu }\nabla p_m\right)+a_{mf}\left(p_f-p_m\right)+a_{mv}\left(p_v-p_m\right);\end{array}$$

(18)

$$\begin{array}{c}\hfill {\varphi }_f{c}_f{\tau }^{{\alpha }_f-1}\frac{{\partial }^{{\alpha }_m}{p}_f}{\partial t^{{\alpha }_m}}=\nabla \left(\frac{k_f}{\mu }\nabla p_f\right)-a_{mf}\left(p_f-p_m\right)+a_{vf}\left(p_v-p_f\right);\end{array}$$

(19)

$$\begin{array}{c}\hfill {\varphi }_v{c}_v{\tau }^{{\alpha }_v-1}\frac{{\partial }^{{\alpha }_m}{p}_v}{\partial t^{{\alpha }_m}}=\nabla \left(\frac{k_v}{\mu }\nabla p_v\right)-a_{mv}\left(p_v-p_m\right)-a_{vf}\left(p_v-p_f\right)\end{array}$$

(20)

where $\tau$ is a coefficient included to maintain dimensional balance in the previous system. Equations (18)–(20) will be referred to as a fractional triple permeability and triple porosity model.

3. Fractional Triple Permeability and Triple Porosity Model

Figure 1a shows how the porous medium is idealized; it is considered that the porous medium is composed of fractures, vuggs and rock matrix where each porous medium interacts with the others. Figure 1b shows how the fluid moves in the porous media and how it reaches the well; as in classical models, it is considered that the matrix leads the fluid to the fractured system and to the vugular system. Besides, the vugular system feeds the fractured system; however, unlike classic models, all three porous media are allowed to feed the well.

Considering the model (18)–(20) and assuming further that: the well penetrates the reservoir completely, the porous medium is single-phase saturated by oil, the interporous flow between each porous medium is pseudo-stable, the density and viscosity are constant, the permeability in each porous medium is constant, the initial pressure $P_i$ is constant throughout the entire reservoir; assuming radial symmetry and that the fluid has an anomalous behavior along each porous medium, that is, including fractional time derivative of Caputo type. The dimensionless model becomes:

$$\begin{array}{c}\hfill {\omega }_m{\tau }_D^{{\alpha }_m-1}\frac{{\partial }^{{\alpha }_m}{p}_{Dm}}{\partial t_D^{{\alpha }_m}}={\kappa }_m\frac{1}{r_D}\frac{\partial }{\partial r_D}\left(r_D\frac{\partial p_{Dm}}{\partial r_D}\right)+{\lambda }_{mf}\left(p_{Df}-p_{Dm}\right)+{\lambda }_{mv}\left(p_{Dv}-p_{Dm}\right),\end{array}$$

(21)

$$\begin{array}{c}\hfill {\omega }_f{\tau }_D^{{\alpha }_f-1}\frac{{\partial }^{{\alpha }_f}{p}_{Df}}{\partial t_D^{{\alpha }_f}}={\kappa }_f\frac{1}{r_D}\frac{\partial }{\partial r_D}\left(r_D\frac{\partial p_{Df}}{\partial r_D}\right)-{\lambda }_{mf}\left(p_{Df}-p_{Dm}\right)+{\lambda }_{vf}\left(p_{Dv}-p_{Df}\right),\end{array}$$

(22)

$$\begin{array}{c}\hfill {\omega }_v{\tau }_D^{{\alpha }_v-1}\frac{{\partial }^{{\alpha }_v}{p}_{Dv}}{\partial t_D^{{\alpha }_v}}={\kappa }_v\frac{1}{r_D}\frac{\partial }{\partial r_D}\left(r_D\frac{\partial p_{Dv}}{\partial r_D}\right)-{\lambda }_{mv}\left(p_{Dv}-p_{Dm}\right)-{\lambda }_{vf}\left(p_{Dv}-p_{Df}\right)\end{array}$$

(23)

where ${\omega }_m=1-{\omega }_f-{\omega }_v$ and ${\kappa }_m=1-{\kappa }_f-{\kappa }_v$, the dimensionless variables are given by

$$\begin{array}{c}\hfill p_{Dj}=\frac{2\pi h(k_m+k_f+k_v)(P_i-p_j)}{Q_0{B}_0\mu },j=m,f,v;\\ \hfill {\omega }_j=\frac{{\varphi }_j{c}_j}{{\varphi }_m{c}_m+{\varphi }_f{c}_f+{\varphi }_v{c}_v},j=m,f,v;\\ \hfill {\kappa }_j=\frac{k_j}{k_m+k_f+k_v},j=m,f,v;\\ \hfill {\lambda }_{ij}=\frac{a_{ij}\mu r_w^2}{k_m+k_f+k_v},i,j=m,f,v;\\ \hfill t_D=\frac{t(k_m+k_f+k_v)}{\mu r_w^2({\varphi }_m{c}_m+{\varphi }_f{c}_f+{\varphi }_v{c}_v)};\\ \hfill r_D=\frac{r}{r_w};\\ \hfill {\tau }_D=\frac{\tau (k_m+k_f+k_v)}{\mu r_w^2({\varphi }_m{c}_m+{\varphi }_f{c}_f+{\varphi }_v{c}_v)};\end{array}$$

the meaning of each variable is given in Abbreviations; $\frac{{\partial }^{{\alpha }_j}{p}_{Dj}}{\partial t_D^{{\alpha }_j}}$ with $j=m,f,v$ is the Caputo fractional derivative defined as

$$\frac{{\partial }^{{\alpha }_j}{p}_{Dj}}{\partial t_D^{{\alpha }_j}}\left(t\right)=\frac{1}{\Gamma (n-{\alpha }_j)}{\int }_0^t\frac{p_{Dj}^{\left(n\right)}\left(\overline{t}\right)d\overline{t}}{{(t-\overline{t})}^{{\alpha }_j-n+1}},$$

(24)

where $p_{Dj}^{\left(n\right)}$ is the usual nth derivative and $n=\left[{\alpha }_j\right]+1$ for ${\alpha }_j\notin \mathbb{N}$, $n={\alpha }_j$ for $\alpha \in \mathbb{N}$ and $\Gamma (·)$ is Euler’s gamma funciton; see [30,31,32].

In general, when the order of the fractional derivative is <1, it is said that the anomalous flow is subdiffusive, while if the order of the fractional derivative is >1, it is said that the anomalous flow is superdiffusive.

By modifying Darcy’s law to introduce the memory phenomenon, Caputo [19] modifies the definition of permeability to maintain dimensional balance, although it introduces the problem of how to measure what it is called “fractional permeability”. In this work, the variable $\tau$ is introduced, with units of speed, independent of the physical parameters of viscosity $\mu$ and permeability $\kappa$; this constant maintains dimensional balance. The dimensional balance variable $\tau$ is kept in the dimensionless model, Equations (21)–(23), through the variable ${\tau }_D$. Section 4 will show the effect of this variable.

From considering the oil field as a porous medium composed of the union of three independent porous media (fractures, vuggs and matrix), added to the assumption that the fluid has an anomalous flow in each porous medium. The proposed model considers the possibility that this anomalous behavior may become different in each porous medium, expressing itself through a different order in the temporal fractional derivative.

In the triple permeability and triple fractional porosity model, constant flow boundary conditions are considered, that is:

$$\begin{array}{c}\hfill \underset{rD\to 1}{lim}{r}_D\left({\kappa }_m\frac{\partial p_{Dm}}{\partial r_D}+{\kappa }_f\frac{\partial p_{Df}}{\partial r_D}+{\kappa }_v\frac{\partial p_{Dv}}{\partial r_D}\right)=-1;\end{array}$$

(25)

$$\begin{array}{c}\hfill p_{Dw}\left(t_D\right)={\left.p_{Dm}(r_D,t_D)\right|}_{r_D=1}={\left.p_{Df}(r_D,t_D)\right|}_{r_D=1}={\left.p_{Dv}(r_D,t_D)\right|}_{r_D=1};\end{array}$$

(26)

$$\underset{r_D\to \infty }{lim}{p}_{Dm}=\underset{r_D\to \infty }{lim}{p}_{Df}=\underset{r_D\to \infty }{lim}{p}_{Dv}=0.$$

(27)

Equation (27) represents the external border condition when considering an infinite reservoir.

Next, the fractional triple permeability and triple porosity model will be analytically solved, Equations (21)–(23), with boundary conditions (25)–(27).

3.1. Laplace Transform

From the classical definition of Laplace transform

$$\overline{p}\left(u\right)=\mathcal{L}\left\{p\left(t\right)\right\}\left(u\right)={\int }_0^{\infty }{e}^{-ut}p\left(t\right)dt,u\ge 0;$$

(28)

where u is the Laplace variable.

It is known that, applying this transformation to the Caputo’s fractional derivative operator, we have

$$\mathcal{L}\left\{\frac{{\partial }^{\alpha }p}{\partial t^{\alpha }}\left(t\right)\right\}=u^{\alpha }\overline{p}\left(u\right)-\sum _{k=0}^{n-1}{u}^{\alpha -k-1}{p}^{\left(k\right)}\left(0\right).$$

(29)

In particular, if $\alpha \in (0,1]$ then $n=1$ and

$$\mathcal{L}\left\{\frac{{\partial }^{\alpha }p}{\partial t^{\alpha }}\left(t\right)\right\}=u^{\alpha }\overline{p}\left(u\right)-u^{\alpha -1}p\left(0\right).$$

(30)

From Equations (29) and (30), it is shown that initial conditions can be incorporated to compute the fractional derivative, which is not the case with other definitions of fractional derivative.

So, the Laplace transform of Equations (21)–(23) are

$$\begin{array}{c}\hfill {\kappa }_m\frac{1}{r_D}\frac{\partial }{\partial r_D}\left(r_D\frac{\partial {\overline{p}}_{Dm}}{\partial r_D}\right)-\left[u^{{\alpha }_m}{\omega }_m{\tau }_D^{{\alpha }_m-1}+{\lambda }_{mf}+{\lambda }_{mv}\right]{\overline{p}}_{Dm}+\\ \hfill {\omega }_m{\tau }_D^{{\alpha }_m-1}\sum _{k=0}^{n-1}{u}^{{\alpha }_m-k-1}{p}_{Dm}^{\left(k\right)}\left(0\right)+{\lambda }_{mf}{\overline{p}}_{Df}+{\lambda }_{mv}{\overline{p}}_{Dv}=0,\end{array}$$

(31)

$$\begin{array}{c}\hfill {\kappa }_f\frac{1}{r_D}\frac{\partial }{\partial r_D}\left(r_D\frac{\partial {\overline{p}}_{Df}}{\partial r_D}\right)-\left[u^{{\alpha }_f}{\omega }_f{\tau }_D^{{\alpha }_f-1}+{\lambda }_{mf}+{\lambda }_{vf}\right]{\overline{p}}_{Df}+\\ \hfill {\omega }_f{\tau }_D^{{\alpha }_f-1}\sum _{k=0}^{n-1}{u}^{{\alpha }_f-k-1}{p}_{Df}^{\left(k\right)}\left(0\right)+{\lambda }_{mf}{\overline{p}}_{Dm}+{\lambda }_{vf}{\overline{p}}_{Dv}=0,\end{array}$$

(32)

$$\begin{array}{c}\hfill {\kappa }_v\frac{1}{r_D}\frac{\partial }{\partial r_D}\left(r_D\frac{\partial {\overline{p}}_{Dv}}{\partial r_D}\right)-\left[u^{{\alpha }_v}{\omega }_v{\tau }_D^{{\alpha }_v-1}+{\lambda }_{mv}+{\lambda }_{vf}\right]{\overline{p}}_{Dv}+\\ \hfill {\omega }_v{\tau }_D^{{\alpha }_v-1}\sum _{k=0}^{n-1}{u}^{{\alpha }_v-k-1}{p}_{Dv}^{\left(k\right)}\left(0\right)+{\lambda }_{mv}{\overline{p}}_{Dm}+{\lambda }_{vf}{\overline{p}}_{Df}=0;\end{array}$$

(33)

where u is the Laplace variable and ${\overline{p}}_{Dj}$ with $j=m,f,v$ is the Laplace transform applied to $p_{Dm}$, $p_{Df}$ and $p_{Dv}$ respectively.

The System (31)–(33) must be treated carefully; if $\alpha \in (0,1]$ then $n=1$ and we would have ${\sum }_{k=0}^{n-1}{u}^{{\alpha }_j-k-1}{p}_{Dj}^{\left(k\right)}(r_D,0)=0$ since $p_{Dj}(r_D,0)=0$ for $j=m,f,v$; however, if ${\alpha }_j\notin (0,1]$ y ${\alpha }_j>0$ with $j=m,f,v$ then additional initial conditions must be included.

3.2. Bessel Function

The solution to each of the equations of the general model (31)–(33) is expressed as a linear combination of modified Bessel functions of first and second type [33]; however, being an infinite site, its solution is expressed as

$$\begin{array}{c}\hfill {\overline{p}}_{Dm}(u,r_D)=A{K}_0\left(\alpha \left(u\right)r_D\right),\end{array}$$

(34)

$$\begin{array}{c}\hfill {\overline{p}}_{Df}(u,r_D)=B{K}_0\left(\alpha \left(u\right)r_D\right),\end{array}$$

(35)

$$\begin{array}{c}\hfill {\overline{p}}_{Dv}(u,r_D)=C{K}_0\left(\alpha \left(u\right)r_D\right).\end{array}$$

(36)

From the substitution of the previous equations, in the model (31)–(33), the matrix system is obtained

$$\left[\begin{array}{ccc}{\kappa }_m{\alpha }^2-m_1& m_2& m_3\\ m_2& {\kappa }_f{\alpha }^2-m_4& m_5\\ m_3& m_5& {\kappa }_v{\alpha }^2-m_6\end{array}\right]\left[\begin{array}{c}A\\ B\\ C\end{array}\right]=-\left[\begin{array}{c}{\omega }_m{\tau }_D^{{\alpha }_m-1}\sum _{k=0}^{n-1}{u}^{{\alpha }_m-k-1}{p}_{Dm}^{\left(k\right)}\left(0\right)\\ {\omega }_f{\tau }_D^{{\alpha }_f-1}\sum _{k=0}^{n-1}{u}^{{\alpha }_f-k-1}{p}_{Df}^{\left(k\right)}\left(0\right)\\ {\omega }_v{\tau }_D^{{\alpha }_v-1}\sum _{k=0}^{n-1}{u}^{{\alpha }_v-k-1}{p}_{Dv}^{\left(k\right)}\left(0\right)\end{array}\right],$$

(37)

where

$$\begin{array}{c}\hfill m_1\left(u\right)=u^{{\alpha }_m}{\tau }_D^{{\alpha }_m-1}{\omega }_m+{\lambda }_{mf}+{\lambda }_{mv},\\ \hfill m_2={\lambda }_{mf},\\ \hfill m_3={\lambda }_{mv},\\ \hfill m_4\left(u\right)=u^{{\alpha }_f}{\tau }_D^{{\alpha }_f-1}{\omega }_f+{\lambda }_{mf}+{\lambda }_{vf},\\ \hfill m_5={\lambda }_{vf},\\ \hfill m_6\left(u\right)=u^{{\alpha }_v}{\tau }_D^{{\alpha }_v-1}{\omega }_v+{\lambda }_{mv}+{\lambda }_{vf}.\end{array}$$

Once again, the order of the fractional derivatives has a big impact, if ${\alpha }_j\in (0,1]$ for $j=m,f,v$ then the matrix Equation (37) turns out to be a homogeneous matrix equation so to obtain non-trivial solutions, it must be satisfied that

$$\begin{array}{c}\hfill {\kappa }_m{\kappa }_f{\kappa }_v{\alpha }^6-\left[{\kappa }_m\left({\kappa }_f{m}_6+{\kappa }_v{m}_4\right)+{\kappa }_f{\kappa }_v{m}_1\right]{\alpha }^4+\\ \hfill +\left[{\kappa }_m\left(m_4{m}_6-m_5^2\right)+\left({\kappa }_f{m}_6+{\kappa }_v{m}_4\right)m_1-{\kappa }_v{m}_2^2-{\kappa }_f{m}_3^2\right]{\alpha }^2\\ \hfill -m_1{m}_4{m}_6+m_1{m}_5^2+m_2^2{m}_6+2m_2{m}_3{m}_5+m_3^2{m}_4=0.\end{array}$$

(38)

Similarly, for $\alpha >0$ and $\alpha \notin (0,1]$ with trivial initial conditions, Equation (38) must be satisfied; for the case with non-trivial initial conditions, the Equation (37) is always trivially solvable if Equation (38) does not hold.

These considerations show the great challenge that arises when approaching a triple permeability and triple porosity model with a fractional temporal derivative; firstly, due to the generalization when considering that the permeability of the matrix is significant and later when adding the fractional temporal derivative in each porous medium, which adds a degree of freedom to the fluid that, as will be seen later, has a significant impact on the well’s pressure.

In the following, it will considered trivial initial conditions; this is $p_{Dj}^{\left(k\right)}=0$ for $j=m,f,v$ and $k=0,1,\cdots ,n-1$.

3.3. Tartaglia–Cardano Equations

Equation (38) is equivalent to a third degree equation, making $\overline{\alpha }={\alpha }^2$. In general, the solution of third degree equation depends of the value of the discriminant $D={(q/2)}^2+{(p/3)}^3$; in particular, for the parameters of a naturally fractured vuggy carbonate reservoir with triple porosity and triple permeability, it follows that $D<0$, which means that roots of (38) are all real numbers. Therefore, applying the Tartaglia–Cardano formulas [34], the solution is

$${\alpha }_k^2=2\sqrt{-\frac{p}{3}}cos\frac{\theta +2k\pi }{3}-\frac{a}{3},cos\theta =\frac{-q/2}{\sqrt{-{(p/3)}^3}},$$

(39)

with $p=b-\frac{a^2}{3}$, $q=\frac{2a^3-9ab+27c}{27}$ and

$$\begin{array}{c}\hfill a=-\frac{\left(1-{\kappa }_f-{\kappa }_v\right)\left({\kappa }_f{m}_6+{\kappa }_v{m}_4\right)+{\kappa }_f{\kappa }_v{m}_1}{\left(1-{\kappa }_f-{\kappa }_v\right){\kappa }_f{\kappa }_v};\end{array}$$

(40)

$$\begin{array}{c}\hfill b=\frac{\left(1-{\kappa }_f-{\kappa }_v\right)\left(m_4{m}_6-m_5^2\right)+\left({\kappa }_f{m}_6+{\kappa }_v{m}_4\right)m_1-{\kappa }_v{m}_2^2-{\kappa }_f{m}_3^2}{\left(1-{\kappa }_f-{\kappa }_v\right){\kappa }_f{\kappa }_v};\end{array}$$

(41)

$$\begin{array}{c}\hfill c=\frac{-m_1{m}_4{m}_6+m_1{m}_5^2+m_2^2{m}_6+2m_2{m}_3{m}_5+m_3^2{m}_4}{\left(1-{\kappa }_f-{\kappa }_v\right){\kappa }_f{\kappa }_v}.\end{array}$$

(42)

3.4. General Solution

Therefore, the solution for pressure comes from replacing (39) in the Equations (34)–(36); consequently, the solution in each porous medium ${\overline{p}}_{Dj}$ for $j=m,f,v$ is rewritten as a linear combination of Bessel function in each solution ${\alpha }_k$, namely

$$\begin{array}{c}\hfill {\overline{p}}_{Dm}(u,r_D)=A_1{K}_0\left({\alpha }_1\left(u\right)r_D\right)+A_2{K}_0\left({\alpha }_2\left(u\right)r_D\right)+A_3{K}_0\left({\alpha }_3\left(u\right)r_D\right),\end{array}$$

(43)

$$\begin{array}{c}\hfill {\overline{p}}_{Df}(u,r_D)=B_1{K}_0\left({\alpha }_1\left(u\right)r_D\right)+B_2{K}_0\left({\alpha }_2\left(u\right)r_D\right)+B_3{K}_0\left({\alpha }_3\left(u\right)r_D\right),\end{array}$$

(44)

$$\begin{array}{c}\hfill {\overline{p}}_{Dv}(u,r_D)=C_1{K}_0\left({\alpha }_1\left(u\right)r_D\right)+C_2{K}_0\left({\alpha }_2\left(u\right)r_D\right)+C_3{K}_0\left({\alpha }_3\left(u\right)r_D\right).\end{array}$$

(45)

By substituting the previous equations in the System (31)–(33), we obtain a matrix system with infinite solutions, similar to the System (37) for ${[A_i,B_i,C_i]}^T$ for $i=1,2,3$. Making $C_1=C_2=C_3=1$ in that matrix system, so it reaches

$$\begin{array}{c}\hfill A_i=\frac{m_3\left({\kappa }_f{\alpha }_i^2-m_4\right)-m_2{m}_5}{m_2^2-\left[\left(1-{\kappa }_f-{\kappa }_v\right){\alpha }_i^2-m_1\right]\left[{\kappa }_f{\alpha }_i^2-m_4\right]},\end{array}$$

(46)

$$\begin{array}{c}\hfill B_i=\frac{-m_3-A_i\left[\left(1-{\kappa }_f-{\kappa }_v\right){\alpha }_i^2-m_1\right]}{m_2};\end{array}$$

(47)

for $i=1,2,3$. Thus, the System (34)–(36) is expressed as follows

$$\begin{array}{c}\hfill {\overline{p}}_{Dm}(u,r_D)=A_1{D}_1{K}_0\left({\alpha }_1{r}_D\right)+A_2{D}_2{K}_0\left({\alpha }_2{r}_D\right)+A_3{D}_3{K}_0\left({\alpha }_3{r}_D\right),\end{array}$$

(48)

$$\begin{array}{c}\hfill {\overline{p}}_{Df}(u,r_D)=B_1{D}_1{K}_0\left({\alpha }_1{r}_D\right)+B_2{D}_2{K}_0\left({\alpha }_2{r}_D\right)+B_3{D}_3{K}_0\left({\alpha }_3{r}_D\right),\end{array}$$

(49)

$$\begin{array}{c}\hfill {\overline{p}}_{Dv}(u,r_D)=D_1{K}_0\left({\alpha }_1{r}_D\right)+D_2{K}_0\left({\alpha }_2{r}_D\right)+D_3{K}_0\left({\alpha }_3{r}_D\right);\end{array}$$

(50)

where $D_1$, $D_2$, $D_3$ are constants determined by the boundary conditions and ${\alpha }_i$ is the positive root of ${\alpha }_i^2$ for $i=1,2,3$.

To obtain the associated constants $D_1$, $D_2$, $D_3$, note that, when replacing the solutions obtained, the Equations (48)–(50) in the boundary conditions, Equations (25) and (26), the following matrix system is reached

$$\left[\begin{array}{ccc}{\alpha }_1{E}_1{K}_1\left({\alpha }_1\right)& {\alpha }_2{E}_2{K}_1\left({\alpha }_2\right)& {\alpha }_3{E}_3{K}_1\left({\alpha }_3\right)\\ (A_1-1)K_0\left({\alpha }_1\right)& (A_2-1)K_0\left({\alpha }_2\right)& (A_3-1)K_0\left({\alpha }_3\right)\\ (B_1-1)K_0\left({\alpha }_1\right)& (B_2-1)K_0\left({\alpha }_2\right)& (B_3-1)K_0\left({\alpha }_3\right)\end{array}\right]\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\end{array}\right]=\left[\begin{array}{c}1/u\\ 0\\ 0\end{array}\right],$$

(51)

where

$$E_i=\left[\left(1-{\kappa }_f-{\kappa }_v\right)A_i+{\kappa }_f{B}_i+{\kappa }_v\right],\mathrm{con}i=1,2,3.$$

(52)

Therefore, the analytical solution in the Laplace space, for pressure at the boundary of the well, $r_D=1$, with dimensionless variables for a single-phase medium with triple permeability and triple porosity is:

$$\begin{array}{c}\hfill {\overline{p}}_{Dw}\left(u\right)=D_1{K}_0\left({\alpha }_1\right)+D_2{K}_0\left({\alpha }_2\right)+D_3{K}_0\left({\alpha }_3\right)\\ \hfill =A_1{D}_1{K}_0\left({\alpha }_1\right)+A_2{D}_2{K}_0\left({\alpha }_2\right)+A_3{D}_3{K}_0\left({\alpha }_3\right)\\ \hfill =B_1{D}_1{K}_0\left({\alpha }_1\right)+B_2{D}_2{K}_0\left({\alpha }_2\right)+B_3{D}_3{K}_0\left({\alpha }_3\right).\end{array}$$

(53)

To calculate pressure in real space, the Stehfest [35] algorithm was used.

Dimensionless pressure considering mechanical skin in the fractured medium, vugular medium and matrix rock, is calculated by

$$p_{Dw}\left(t_D\right)={\left.p_{Dm,Df,Dv}\right|}_{r_D=1}-s_{m,f,v}{\left.r_D\frac{\partial p_{Dm,Df,Dv}}{\partial r_D}\right|}_{r_D=1}.$$

(54)

By using this condition instead of (26), constants $D_1$, $D_2$ and $D_3$ are determined and the following matrix system is reached

$$\left[\begin{array}{ccc}{\alpha }_1{E}_1{K}_1\left({\alpha }_1\right)& {\alpha }_2{E}_2{K}_1\left({\alpha }_2\right)& {\alpha }_3{E}_3{K}_1\left({\alpha }_3\right)\\ \begin{array}{c}(A_1-1)K_0\left({\alpha }_1\right)+\\ {\alpha }_1{K}_1\left({\alpha }_1\right)(s_m{A}_1-s_v)\end{array}& \begin{array}{c}(A_2-1)K_0\left({\alpha }_2\right)+\\ {\alpha }_2{K}_1\left({\alpha }_2\right)(s_m{A}_2-s_v)\end{array}& \begin{array}{c}(A_3-1)K_0\left({\alpha }_3\right)+\\ {\alpha }_3{K}_1\left({\alpha }_3\right)(s_m{A}_3-s_v)\end{array}\\ \begin{array}{c}(B_1-1)K_0\left({\alpha }_1\right)+\\ {\alpha }_1{K}_1\left({\alpha }_1\right)(s_f{B}_1-s_v)\end{array}& \begin{array}{c}(B_2-1)K_0\left({\alpha }_2\right)+\\ {\alpha }_2{K}_1\left({\alpha }_2\right)(s_f{B}_2-s_v)\end{array}& \begin{array}{c}(B_3-1)K_0\left({\alpha }_3\right)+\\ {\alpha }_3{K}_1\left({\alpha }_3\right)(s_f{B}_3-s_v)\end{array}\end{array}\right]\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\end{array}\right]=\left[\begin{array}{c}1/u\\ 0\\ 0\end{array}\right]$$

(55)

giving the solution for the pressure

$${\left.{\overline{p}}_{Dw}\right|}_{C_D=0}=D_1{K}_0\left({\alpha }_1\right)+D_2{K}_0\left({\alpha }_2\right)+D_3{K}_0\left({\alpha }_3\right)+s_v(D_1{\alpha }_1{K}_1\left({\alpha }_1\right)+D_2{\alpha }_2{K}_1\left({\alpha }_2\right)+D_3{\alpha }_3{K}_1\left({\alpha }_3\right)).$$

(56)

To include the effect of storage in wellbore, the following equation is used:

$${\overline{p}}_{Dw}=\frac{1}{u^2{C}_D+1/{\left.{\overline{p}}_{Dw}\right|}_{C_D=0}}.$$

(57)

4. Results and Discussion

This section will present several synthetic results of solving the fractional triple permeability and triple porosity model, Equations (21)–(23), subject to the initial and boundary conditions (25)–(27). The results show the value of the pressure at the boundary of the well in real space $p_{Dw}$ and its respective Bourdet derivative, $p_{Dw}^{\prime }$ [8]; pressure results will be shown with a solid line, while its derivative will be with a dashed line; both results will be drawn with the same color.

Figure 2 shows the results of the triple permeability and triple porosity fractional model compared to other models in the literature. In this figure, a skin value of $s=s_m=s_f=s_v=1.5$ and a well storage coefficient of $C_D=100$ are considered. Figure 2a shows the results considering a non-anomalous flow, i.e., ${\alpha }_m={\alpha }_f={\alpha }_v=1$; It is shown how ${\kappa }_m\to 0$ recovers the triple porosity and double permeability model [10] and; moreover, results of double porosity models [22]. Figure 2b shows the results of considering an anomalous flow; in particular, it is highlighted the similarity between the results of Razminia et al. [22] considering the Euclidean case $\theta =0$ and the results of the triple permeability and triple porosity fractional model considering an anomalous flow throughout the whole reservoir, i.e., $\alpha ={\alpha }_m={\alpha }_f={\alpha }_v$.

The next results will show the effect of considering that the fluid has an anomalous behavior in each porous medium independently and as a whole. The following results consider a skin and null well storage effect; that is, $s_j=0$ for $j=m,f,v$ and $C_D=0$.

Considering the fluid flow in reservoirs with triple porosity and non-anomalous flow; the pressure deficit $p_{Dw}$ can be separated into three main phases: short-time phase, where the well storage effect is shown through the straight line behavior with slope 1, followed by the skin effect and ending with the so-called fractures dominated phase, behavior that is exhibited in short times if the previous effects were null. The transient phase, where each porous medium dominates the pressure deficit in the well, starting with the phase dominated by the fracture system, until it loses its domain and yields it to the vugular system and ends with the phase dominated by the matrix system; in particular, the effect of the interporous flow coefficients indicate the moment in which each porous medium ends its domain, being reflected by a “V”-shaped behavior, forming, as a whole, what is called “W”-shaped behavior when said behaviors do not intersect. Finally, the stability phase, when the stability phase of the matrix system is reached and the pressure deficit grows at a constant rate, being reflected in its derivative through a horizontal line.

4.1. Case Study 1

Figure 3 considers that the fluid has an anomalous behavior only through the fractured system; that is, in Equations (21)–(23), the order of the fractional derivative for the vugular system and the rock matrix is 1. The Figure shows how the effect of the anomalous flow influences the phenomenon from short times, phase dominated by fractures, where for superdiffusive flows, it is shown with greater force, compared to the case of subdiffusive flows, where the effect of the fractional time derivative is barely distinguishable. See Appendix A.

The transient phase is particularly insensitive to anomalous flow from the fractured system; however, an additional phenomenon can be observed between the transient phase and the stability phase for subdiffusive flows. This phenomenon is characterized by a drop in the rate of pressure increase that affects the second crest of the characteristic “W”-shaped behavior when ${\alpha }_f\to 0$ affecting the stability phase, both in its appearance and in the growth rate.

Throughout the evolution of the phenomenon, the fractional time derivative impacts the effect of interporous flows; in particular for the anomalous flow in the fractured system, the interporous flow coefficient between the matrix and the fracture, ${\lambda }_{mf}$, is affected. Therefore, for subdiffusive flows, the feeding to the fractured system is slower, which explains the pressure drop instead of reaching stability.

Figure 4 considers that the fluid has an anomalous behavior only through the vugular system, that is, in Equations (21)–(23), the order of the fractional derivative for the fracture system and the rock matrix is 1. The figure shows the impact of the anomalous vugular flow beginning in the phase dominated by fractures whose behavior is analogous to the case of anomalous flow in the fractured system, see Appendix A; however, its effect is longer lasting, affecting the transient phase.

During the transient phase, the anomalous flow influences the characteristic “W”-shaped behavior; in particular, it impacts on the effect of the interporous flow coefficient between the matrix and the vugular system, ${\lambda }_{mv}$, increasing its effect for superdiffusive flows and inhibiting it for subdiffusive cases; which is shown in the first crest of the “W” shaped behavior.

The final part of the transient phase is immune to the anomalous flow of the vugular system; however, and analogously to the case of anomalous flow in the fractured system, a drop in the rate of increase in pressure is observed that affects the second crest of the characteristic “W”-shaped behavior in the case of subdiffusive flows affecting the stability phase.

Figure 5 shows the case of an anomalous flow in the matrix, making the order of the fractional time derivative for the fractured system and the vugular system 1. As expected from the asymptotic behavior in short times, see Appendix A, the effect of the anomalous flow of the matrix is greater for the phase dominated by fractures compared to the effect shown in the previous cases.

On the other hand, the effect that the anomalous flow of the matrix has on the interporous flow coefficients, ${\lambda }_{mf}$ and ${\lambda }_{mv}$, considerably alters the characteristic “W”-shaped behavior, completely altering the transient phase; this affects the shape of the crests of the “W” -shaped behavior, as well as the time in which they occur.

In particular, for superdiffusive flows, the stability phase is delayed due to the overexcitation received by the porous media coming from the matrix; while for subdiffusive flows, the same phenomenon of drop in the pressure increase rate seen in previous cases is observed.

Figure 6 shows the results of considering an anomalous flow in the whole porous medium composed of the fractured, vugular and matrix systems; in particular, the anomalous behavior is considered to be congruent in each porous medium; that is, $\alpha ={\alpha }_f={\alpha }_v={\alpha }_m$ in (21)–(23) where the combined effect of the fractional time derivative of each porous medium will be seen. Like in the case of anomalous flow in each porous medium, it affects the phase dominated by fractures; when considering the anomalous flow in the whole porous medium, it will show, with greater impact, its effect in short times, see Appendix A.

As might be expected, from the effect caused by the anomalous vugular flow and increased by the anomalous flows from the other porous media, the phase dominated by fractures culminates with a greater impact that has repercussions in the beginning of the transient phase.

The transient phase has a combined impact by the anomalous flow of each porous medium, resulting from the impact to the interporous flow coefficients ${\lambda }_{ij}$ with $i,j=f,v,m$.

Finally, the stability phase is reached promptly for superdiffusive casesm because the combined effect of anomalous flows generates an increase in the rate of increase pressure in contrast to the effect of independent anomalous flows where the rate of increase in the pressure tends to be the same as in the non-anomalous case; while for subdiffusive cases, the stability phase is delayed, showing a behavior as in the cases already described.

4.2. Case Study 2

The previous case study showed the effect of an anomalous flow isolated to each porous medium and an anomalous flow in the entire porous medium, combined with a remarkably significant permeability in the matrix. In this case study, the same phenomenon will be described, highlighting the effect of permeability in the matrix, contrasting the previous results with the results of an almost negligible permeability in the matrix, that is, ${\kappa }_m\approx 0$.

From the equations in Appendix A, an almost negligible permeability in the matrix strongly affects the fracture-dominated phase by pronouncing the pressure velocity in the wellbore; furthermore, the duration of this phase is also extended, as can be seen in Figure 7, Figure 8, Figure 9 and Figure 10.

Figure 7 shows the combined effect of almost negligible permeability with the effect of anomalous flow in the fractured system. In particular, having a nearly impermeable matrix increases the effect generated by the interporous flow coefficients, namely ${\lambda }_{mf}$ and ${\lambda }_{mv}$, where the characteristic “W”-shaped behavior is more remarkable.

By having a more marked “W”-shaped behavior, the effect created by the anomalous flow from the fractured system has on the effects generated by the interporous flow coefficients is even more remarkable showing, for subdiffusive flows, the phenomenon described in the case of study 1; whereas, for superdiffusive flows, a trough mitigation typical of the interporous flow effect ${\lambda }_{mf}$ is shown when ${\alpha }_f\to 2$, precipitating the stability phase.

Figure 8 shows more strongly the effect created by an anomalous flow in the fractured system; where the phase dominated by fractures, being more marked, extends its duration by increasing the first crest that makes up the “W” -shaped behavior characteristic of the transitory phase, which directly impacts the trough that follows; however, the anomalous flow does not affect all the “W”-shaped behavior, so the end of the transient phase and the stability phase is as in the previous case.

As in Figure 5, Figure 9 shows the effect of anomalous flow in the matrix. In particular, when the effects of permeability in the matrix in non-anomalous cases are compared; it is observed that by decreasing the permeability value in the matrix, in addition to impacting the phase dominated by fractures, the effect on interporous flows increases, pronouncing the characteristic “W”-shaped behavior in the transient phase; consequently, the time in which the stability phase is reached is delayed.

Thus, the anomalous flow in the matrix enhances the effects generated by the permeability of the matrix by altering the way in which the fractured system and the vugular system are fed by the matrix. In particular, superdiffusive flows accelerate at the moment when the crests and troughs of the characteristic “W”-shaped behavior shows up by overstimulating said crests and troughs; while for subdiffusive behaviors, the flow follows the same behavior described in the previous case.

Figure 10 shows the equivalent anomalous flow effect in each porous medium, where, as discussed in the previous case, a combined effect is generated from considering an anomalous flow in each porous medium; besides, considering an almost negligible permeability in the matrix, shows an effect that softens the phase dominated by fractures, impacting, as a consequence, the time in which the transient phase begins.

In the same way as has been discussed above, the variation in permeability influences the effect caused by the interporous flow coefficients, so when adding the anomalous flow phenomenon, the crests and troughs that make up the characteristic behavior are overstimulated.

Finally, the stability phase is invariable to the effect on the variation of the permeability in the matrix, so the anomalous flow effect is analogous to that discussed in the previous case.

Figure 11 shows the results of assuming an anomalous flow in the whole reservoir, considering a non-trivial dimensionless coefficient ${\tau }_D\ne 1$. As can be confirmed in the appendix, the dimensionless coefficient ${\tau }_D$ strongly impacts the phase dominated by fractures and consequently the subsequent phases, appearing mainly as a delay.

As observed comparing the previous case studies, by eliminating the assumption of negligible matrix permeability and considering it significant, ${\kappa }_m\ne 0$, a degree of freedom is added with respect to the existing models in the literature, allowing modeling a broader set of heterogeneous reservoirs than those that can be obtained with triple porosity and double permeability models, making ${\kappa }_m=0$; triple porosity and single permeability models, making ${\kappa }_f=0$ for primary flow in the fractured system and ${\kappa }_f=1$ for primary flow in the vugular system; and double porosity models. On the other hand, each case study explores anomalous flows not limited to subdiffusive flows as in the existing literature; additionally, the order of the fractional derivative is allowed to be different in each porous medium, showing particular results of each porous medium, exploiting the idea that, when considering a porous medium with anomalous flow, each porous medium may have particular anomalous flows and different from those of the surrounding porous media.

5. Conclusions

A new model is proposed to capture the high heterogeneity in naturally fractured vuggy reservoirs saturated with one fluid, oil, which relaxes the usual assumptions in classical models by considering that the effect of permeability in the matrix is significant; in addition, all media porous that make up the reservoir (fractured system, vugular system and matrix) contribute to the well feeding. Furthermore, the phenomenon of anomalous flows is considered, incorporating the Caputo type temporary fractional derivative with orders in $(0,2)$, encompassing the phenomena known as subdiffusive and superdiffusive; constants are added to the resulting system, in order to maintain dimensional balance and ensure that the physical conclusions obtained are consistent. The model is solved analytically considering an infinite reservoir in Laplace space, solving an equivalent third degree equation with the Tartaglia–Cardano formulas, and obtaining results in real space, applying the Stefhest method. The results showed the anomalous flow effect that each porous medium provides separately and as a whole to the phenomenon of transient pressure. The impact that each anomalous flow has in the phase dominated by fractures is exposed, the effect that the anomalous flow of each porous medium has on each part of the transitory phase is highlighted, reaching to pronounce or soften each crest and trough in the characteristic “W”-shaped behavior, depending on the porous medium that has the anomalous behavior. Likewise, the phenomenon of decay in fluid velocity generated by subdiffusive flows is shown, creating a stability phase markedly different from the classical case, while superdiffusive flows show a stability phase similar to the non-anomalous case when the anomalous phenomenon is restricted to a single porous medium, while when the anomalous phenomenon is shown in the complete reservoir, the stability phase is altered, showing a higher velocity in the fluid pressure.