Spatial Fractional Darcy’s Law on the Diffusion Equation with a Fractional Time Derivative in Single-Porosity Naturally Fractured Reservoirs

: Due to the complexity imposed by all the attributes of the fracture network of many naturally fractured reservoirs, it has been observed that ﬂuid ﬂow does not necessarily represent a normal diffusion, i.e., Darcy’s law. Thus, to capture the sub-diffusion process, various tools have been implemented, from fractal geometry to characterize the structure of the porous medium to fractional calculus to include the memory effect in the ﬂuid ﬂow. Considering inﬁnite naturally fractured reservoirs (Type I system of Nelson), a spatial fractional Darcy’s law is proposed, where the spatial derivative is replaced by the Weyl fractional derivative, and the resulting ﬂow model also considers Caputo’s fractional derivative in time. The proposed model maintains its dimensional balance and is solved numerically. The results of analyzing the effect of the spatial fractional Darcy’s law on the pressure drop and its Bourdet derivative are shown, proving that two deﬁnitions of fractional derivatives are compatible. Finally, the results of the proposed model are compared with models that consider fractal geometry showing a good agreement. It is shown that modiﬁed Darcy’s law, which considers the dependency of the ﬂuid ﬂow path, includes the intrinsic geometry of the porous medium, thus recovering the heterogeneity at the phenomenological level.


Introduction
The modeling of physical phenomena is always a complex task in which hypotheses and idealizations are used in order to implement already studied laws or similar models.In particular, modeling the transient pressure behavior in naturally fractured reservoirs is a task that has been studied and worked on by various researchers over decades, where the heterogeneity of the porous medium and anomalous fluid flow has been a challenge for which more complex models have been generated and implemented, obtaining results that constantly improve the understanding of the phenomenon.
Considering the complexity of the porous medium, Chang and Yortsos [1] are the pioneers that considered the porous medium as a fractal system, allowing them to describe reservoirs with spatial disorder and, therefore, a complex fluid flow path.
In recent years, fractional calculus has become a useful tool that, applied to diffusiontype problems, explains the behavior of anomalous flow (where the movement of the fluid does not have a Brownian-type behavior) by demonstrating that the fluid has memory [2].
The study of Metzler et al. [3] has become an important reference when studying the phenomenon of anomalous diffusion in a medium with a fractal structure, allowing endless researchers to explore the implications considering different approaches [4][5][6][7].
Although all these models have allowed a more precise understanding of the phenomenon in non-homogeneous reservoirs; the task of understanding this complex phenomenon has also included modifications to Darcy's law.When considering non-Darcy flows, i.e., a sub-diffusive process, to include the memory effect, the fractional derivative approach is used.
To include the memory effect in Darcy's law, Caputo's fractional time derivative is used in Darcy's law modeling, to account for the effect of a decrease in permeability with time [8].Raghavan [9] shows that using the properties of the fractional derivative, one can translate the fractional time derivative into the continuity equation.
Several variants of the temporal fractional derivative have been applied to Darcy's law.Chang et al. [10] replace the spatial derivative in Darcy's law with the Riemann-Liouville fractional derivative to quantify the spatial path dependence of a fluid flow.El-Amin [11] constructs the fractional mass equation and combines it with variants of Darcy's law that include adding the fractional time derivative and replacing the space derivative with Caputo's fractional derivative.Chang and Sun [12] replace the time derivative with Caputo's fractional time derivative to describe the heavy tail decay at long times and the space derivative with the Riemann-Liouville fractional time derivative to describe the dependence of non-local concentration change on a wide range of spatial areas.
Likewise, Caputo and Plastino [13], considering infinite media, add a term proportional to the spatial fractional derivative of the Caputo type to the classic Darcy's law, which represents the effect of spatial memory, i.e., the pressure gradient investigated from the initial point to the measured point.
Obembe et al. [14], in their excellent review work, show a model of Darcy's law, among others, where the temporal fractional derivative is applied to a complex combination of a spatial fractional derivative and complementary temporal fractional derivative; that reflects the presence of flow hindrances, while the spatial fractional derivatives consider flow buffers.
Finally, in order to consider the intrinsic geometry of the porous medium, Cloot and Botha [15] replace the spatial derivative in Darcy's law with Weyl's fractional derivative; they consider that the fluid flow at a given point is governed not only by the properties of the pressure field at a specific position but also depends on the global spatial distribution of that field.
Therefore, in order to model the transient pressure behavior in a naturally fractured reservoir, including spatial and temporal memory, a new model will be proposed that incorporates both the fractional time derivative in the Caputo sense, in the continuity equation, and the fractional Weyl derivative, in Darcy's law.This work is organized as follows: Section 2 describes the essential mathematical tools used; Section 3 shows the development of the model, explaining the modified Darcy's law used and the spatial fractional diffusion model developed; Section 4 shows the results of solving this model and the effect of every parameter of the model, Section 5 compares the proposed model with models that consider the porous medium as fractal; in the end, Section 6 describes the main conclusions reached.

Mathematical Background
This section presents the mathematical tools that will be applied throughout this work.For more details, see [16][17][18].Definition 1.Let t 0 < ∞.The Riemann-Liouville fractional integral RL I α t 0 + y of order α ∈ R is defined by where Γ(•) is Euler's gamma function.
The Caputo fractional derivative, expressed from the fractional integral, is defined as follows: Definition 2. Let t 0 < ∞.For α ∈ R, the Caputo fractional derivative, C D α t 0 + y is defined by where D = d/dt and n ∈ N with n − 1 < α ≤ n.
Proposition 1.The Caputo fractional derivative provides an inverse operator to the Riemann-Liouvile fractional integral, that is The above definitions consider t 0 < ∞; however, similar definitions exist on the whole axis R.
Indeed, the Weyl fractional integral is interesting; however, care should be taken when using this definition because it may not be applied to all functions.
Partcularly, if y(t) is integrable on any finite subinterval of J = [0, ∞), and, if y(t) behaves like t −µ for t large, then the Weyl fractional integral of y of order β will exist if 0 < β < µ.Definition 4. If y(t) is a function for which W −β y(t) exists and has n continous derivatives; then, the Weyl fractional derivative of y of order ν ∈ R is defined by where D = d/dt and n ∈ N with n − 1 < β ≤ n.
Before finishing this long list of definitions, some interesting properties of fractional operators will be shown.
Both, the fractional Weyl integral and the fractional Weyl derivative satisfy the following:

Model Development
In this section, the proposed Darcy's law will be presented, which integrates the fractional derivative of Weyl, and the flow equation with spatial fractional Darcy flow will be constructed.
We consider a fully penetrated well in an infinite porous media with a single porosity system, i.e., a Type I fractured system of Nelson [19], constant initial pressure, permeability, density and viscosity.Furthermore, we consider that the fluid flow at a given point is governed not only by the properties of the pressure field at that specific position but also depends on the global spatial distribution of that field, and lastly, we consider that the radial symmetry is valid.The model in dimensionless variables with a fractional time derivative for the fluid transfer equation, resulting from the combination of the continuity equation and Darcy's law, in radial coordinates and considering a Euclidean porous medium, is given by whereas, when considering a fractal porous medium, that is, a fractal reservoir, the model is [20].In Equations ( 7) and (8), the fractional time derivative is the Caputo fractional derivative, For both models, the following initial and boundary conditions are considered.
• Inner boundary condition lim • Outer boundary condition lim The dimensionless variables are defined as follows: where τ is a constant introduced for the purpose of maintaining dimensional balance.
The proposed model considers a slightly compressible liquid; that is, all the complexities due to the multiphase gas-oil interactions are not taken into consideration [21].Furthermore, the model presented in this paper does not consider petrophysical properties dependent on the stress state, which even for a Nelson Type I reservoir, imposes a challenge [22].

Spatial Fractional Darcy's Law
Figure 1 shows the solution of Equations ( 7) and ( 9)- (11), i.e., the behavior of the pressure drop, p D , throughout space for different times and values of the fractional time derivative order, α. Figure 1 shows the decreasing behavior of the pressure drop in space and, therefore, shows that it is valid to apply the Weyl fractional derivative to p D .
For that reason, it is proposed to modify Darcy's law by substituting the spatial derivative for the Weyl fractional derivative, obtaining the following spatial fractional Darcy's law: p D and, therefore, the classical Darcy's law is recovered.

Spatial Fractional Diffusion Model
Given the continuity equation, the first equation in (7), the traditional Darcy's law, q D , is replaced by the spatial fractional Darcy's law, Equation (12), where the spatial fractional diffusion model is obtained, namely As a consequence of the spatial fractional Darcy's law, the inner boundary condition is also modified, namely lim Therefore, the spatial fractional diffusion model to be solved is the one constituted by Equations ( 9), ( 11), ( 13) and (14).
With these modifications, it is intended that the fluid flow at a given point in the porous medium be governed by the global spatial distribution of the pressure field and not only by the behavior in the direct neighborhood of the pressure around that point.

Results
In this section, the solution of the spatial fractional diffusion model will be shown, as well as its behavior when its parameters vary, and such behavior will be compared with the model that considers a fractal reservoir.
The spatial fractional diffusion model was solved using the finite difference method, considering an implicit scheme in time and Crank-Nicholson in space.
Figure 2 shows the numerical solution of solving the spatial fractional diffusion model considering the classical time derivative, α = 1, and varying the order in the fractional Darcy's law, β.The solution in the well, p wD = p D (r D = 1, t D ), and its Bourdet semilogarithmic derivative, p wD [23], are shown.In Figure 2, two important features can be described when modifying β in the spatial fractional Darcy's law.
At short times, the pressure drop is lower than the pressure drop of the classical case, β = 1, which corresponds to a higher order in the spatial fractional Darcy's law; while if β is lower than one (i.e., classical case), a higher pressure drop will be obtained.However, in the Bourdet derivative, β in the spatial fractional Darcy's law does not seem to have a remarkable impact, except by keeping a power-law behavior with a slope of 0.5.
At long times, when β < 1, the pressure drop is lower than that in the classic case and its respective Bourdet derivative shows a decrease following a power-law behaviour, showing a greater connectivity of pores that provide preferential flow paths; while for β > 1, the pressure drop is greater than in the classic case following a power-law behaviour, and the same happens with the pressure derivative.Both power-law behaviors are parallel, evidencing the creation of flow buffers in the porous medium.
Figure 3 shows the combined effect of the fractional temporal derivative and the spatial fractional Darcy's law, that is, for different values of β in the space fractional Darcy's law, on the pressure drop and its Bourdet derivative for different values of the order of the fractional time derivative, α.It is necessary to highlight that the results shown in Figure 3 demonstrate the compatibility of two different definitions of fractional derivative; that is, Caputo's fractional derivative for the temporal derivative and Weyl's fractional derivative for the spatial derivative in Darcy's law, which allows better modeling of a complex phenomenon by making each fractional derivative capture different characteristics of this phenomenon.
Figure 3a shows for β < 1 the effect of varying the order of the fractional time derivative, α.In the case with the integer time derivative, α = 1, the effect of spatial fractional Darcy's law creates a larger short-time pressure drop that quickly loses its effect, making the long-time pressure drop smaller the smaller α is; this is a consequence of a greater connectivity of pores that provide preferential flow paths.In the case of β = 0.5, Figure 3a, we can observe that the influence of α is not as big as that observed in the pressure derivative, where we observe power-law behaviors at late times, which could be an indication of flow restrictions.
Figure 3b shows for β > 1 the effect of varying the order of the fractional time derivative, α, on the pressure deficit and its Bourdet derivative.For the case with an integer time derivative, α = 1; at short times, the spatial fractional Darcy's law causes the fluid flow velocity to have a behavior similar to the case with the classical Darcy's law, i.e., a linear flow behavior; however, the pressure drop is smaller.Subsequently, for long times, the pressure drop is lower than in the classical case, the lower α is; this is caused by the displacement of dead pores that creates preferential flow paths, which is reflected in the increase in flow velocity and, therefore, in the decrease of the pressure drop.By incorporating the memory effect of the fluid flow, the dependence of the previous pressure on the subsequent pressures is observed; that is, the flow of the fluid receives an increase in the phenomenon that is reflected in the decrease of the pressure drop and its Bourdet derivative.
It is important to point out that in all cases with β = 1.5 and α less or equal to one, Figure 3b, at late times, there is parallel power-law behavior for both the pressure drop and semilog pressure derivative.This behavior is similar to the fractal behavior, and both cases are expressions of anomalous diffusion; in this case, there is a super-diffusion because of the memory effect, and in the fractal case, there is sub-diffusion.
Figure 4 shows the effect of δ β,D , the dimensionless variable introduced in the spatial fractional Darcy's law to maintain the dimensional balance in the fractional flow equation, for the cases β < 1 and β > 1 with α = 1.The effect of δ β,D is remarkable throughout the whole phenomenon.For β < 1 and in short times, it is observed that the greater the value of δ β,D , the greater the value of p wD , also affecting the transition phase, presenting itself at different times for each value of δ β,D and finally reaching the same value of p wD in the stability phase.Bourdet's derivative shows that although in short times, the greater the value of δ β,D the greater p wD will be and the faster the decrement in p wD to such a degree that upon reaching the stability phase, it will be shown that the higher the value of δ β,D , the lower the value of p wD , also showing a parallel behavior for all the values of δ β,D .
On the other hand, for β > 1 at short times, it is observed that in the same way as in the previous case, the greater δ β,D , the greater the value of p wD , changing the growth rate in the transient phase and subsequently having a constant growth rate for each value of δ β,D .This behavior is maintained in the Bourdet derivative, where, in short times, the higher the value for p wD , the higher the δ β,D , followed by a parallel growth for the different values of δ β,D , changing the growth rate during the transient phase until reaching a constant and parallel growth for the different values of δ β,D .

Discussion
In this section, the spatial fractional diffusion model, derived from the equation of continuity with a fractional time derivative and the proposed spatial fractional Darcy's law, and the results from the fractal model are compared.
The purpose is to compare, through the results shown, the approaches considered by both models.The spatial fractional diffusion model shows a great agreement with the data of the model that considers a fractal reservoir, both for the pressure drop and its Bourdet derivative; thereby, it can be considered that by including the spatial fractional Darcy's law, the heterogeneity of the porous medium is recovered through the use of the fractional approach.Furthermore, the above makes sense considering that the spatial fractional Darcy's law takes into account the geometry of the porous medium in two ways; on the one hand, the geometry of the porous medium is implicitly included through the parameter in the Weyl fractional derivative, β; while, on the other hand, it is explicitly included through the δ parameter, which also helps with dimensional balance.
With the above, it is shown that the spatial fractional Darcy's law, by including the dependence on the fluid flow path, i.e., spatial memory, also takes into account the intrinsic geometry of the porous medium.
The proposed model can be part of a robust deep learning model, such as the one suggested by [24], for the automated analysis of pressure tests.

Conclusions
Considering the hypotheses that the Type I naturally fractured reservoir of Nelson is embedded in an infinite porous medium, with a slightly compressible fluid and stressindependent petrophysical properties, we have the following conclusions: 1.
Fluid flow at a given point is governed not only by the properties of the pressure field at that specific position but also depends on the global spatial distribution of that field; it is proposed to modify Darcy's law by replacing the spatial derivative with the Weyl fractional derivative.

2.
The proposed model, with the Caputo fractional time derivative and Weyl fractional derivative on the gradient, proves that two different definitions of the fractional derivative can be compatible.

3.
The results show that the spatial fractional Darcy's law can, on the one hand, reflect the preferential flow paths, and on the other hand, the creation of flow buffers during the phenomenon, whereas the time-fractional derivative incorporates the memory effect of the fluid flow.4.
The proposed model not only resembles the results obtained in models that incorporate fractal geometry, but it also incorporates the intrinsic geometry of the porous medium and, therefore, allows recovering the heterogeneity of the porous medium.This is important because it is confirmed for the first time that both the fractal and fractional approaches represent two possible ways of capturing anomalous diffusion.

Nomenclature
The following abbreviations are used in this manuscript:

Figure 1 .
Figure 1.Spatial pressure behavior at different times, considering the Euclidean porous medium.
= −W β p D and δ β,D = r β < 2 and δ is a constant term included to maintain dimensional balance in the flow equation.Note that setting β = 1,

Figure 2 .
Figure 2. Effect of spatial fractional Darcy's law with Weyl's fractional derivative, considering the classical time derivative α = 1 with δ β,D = 1.The solid lines denote the pressure drop in the well, p wD , and the dashed lines their Bourdet derivative, p wD .

Figure 3 .
Figure 3.Effect on the pressure drop and its Bourdet derivative of the combined effect of the fractional temporal derivative and the spatial fractional Darcy's law.(a) For β < 1, the behavior when varying the order of the fractional time derivative.(b) For β > 1, the behavior when varying the order of the fractional time derivative.In both cases, δ β,D = 1.

Figure 4 .
Figure 4. Effect of δ β,D in pressure drop and Bourdet derivative of the spatial fractional diffusion model with α = 1.0.
p D Dimensionless pressure of the porous medium p wD Dimensionless pressure at the well t D Dimensionless time r D Dimensionless distance q D Dimensionless flow rate of fluid of the porous medium α Order of the fractional time derivative τ D Auxiliary constant to get dimensional balance in the dimensionless fractional differential equation d m f Mass fractal dimension θ Conductivity index P i Constant initial pressure (Pa) p Pressure in the porous media (Pa) φ Porosity of the porous media (m 3 /m 3 ) c Compressibility of the porous medium (Pa −1 ) κ Permeability of the porous medium (m 2 ) µ Dynamic viscosity of the fluid (Pa•s) t Time (s) r Distance variable (m) r w Well radius (m) h Reservoir height (m 2 ) Q 0 Flow rate of extracted fluid (m 3 •s −1 ) B 0 Oil formation volume factor, RB/STB q Flow rate of fluid per unit area of the porous medium (m/s) τ Auxiliary constant to get dimensional balance in the fractional differential equation (s) q β,D Dimensionless spatial fractional Darcy's law δ β,D Auxiliary constant to get dimensional balance in the dimensionless spatial fractional Darcy's law δ Auxiliary constant to get dimensional balance in the spatial fractional Darcy's law (m)