Modeling of 1D Anomalous Diffusion in Fractured Nanoporous Media
Modélisation de la diffusion anormale en 1D dans des milieux nanoporeux fracturés
Colorado School of Mines, Petroleum Engineering Department, 1600 Arapahoe Street, Golden, CO
80401  USA
email: aalbinal@mymail.mines.edu  rholy@mymail.mines.edu  hsarak@mymail.mines.edu  eozkan@mymail.mines.edu
^{*} Corresponding author
Received:
31
May
2015
Accepted:
25
May
2016
Revised:
25
November
2015
Fractured nanoporous reservoirs include multiscale and discontinuous fractures coupled with a complex nanoporous matrix. Such systems cannot be described by the conventional dualporosity (or multiporosity) idealizations due to the presence of different flow mechanisms at multiple scales. More detailed modeling approaches, such as Discrete Fracture Network (DFN) models, similarly suffer from the extensive data requirements dictated by the intricacy of the flow scales, which eventually deter the utility of these models. This paper discusses the utility and construction of 1D analytical and numerical anomalous diffusion models for heterogeneous, nanoporous media, which is commonly encountered in oil and gas production from tight, unconventional reservoirs with fractured horizontal wells. A fractional form of Darcy’s law, which incorporates the nonlocal and hereditary nature of flow, is coupled with the classical mass conservation equation to derive a fractional diffusion equation in space and time. Results show excellent agreement with established solutions under asymptotic conditions and are consistent with the physical intuitions.
Résumé
Des réservoirs nanoporeux fracturés comprennent des fractures à échelles multiples et discontinues couplées à une matrice nanoporeuse complexe. Ces systèmes ne peuvent pas être décrits par les modélisations conventionnelles à double porosité (ou multiporosité) du fait de la présence de différents mécanismes de flux à multiples échelles. Des approches de modélisation plus détaillées, telles que les modèles de Réseau de Fracture Discret (DFN, Discrete Fracture Network) souffrent de même des exigences extensives en matière de données, dictées par la complexité des échelles de flux, qui nuisent éventuellement à l’utilité de ces modèles. Le présent article traite de l’utilité et de la construction de modèles analytiques et numériques de diffusion anormaux en 1D pour des milieux hétérogènes nanoporeux, qui se rencontrent communément dans la production de pétrole et de gaz, pour les réservoirs étanches, non conventionnels avec des puits horizontaux fracturés. Une forme fractionnée de la loi de Darcy, qui comprend la nature nonlocale et héréditaire du flux, est reliée à l’équation de conservation de la masse classique afin d’obtenir une équation de diffusion fractionnée dans l’espace et le temps. Les résultats montrent une parfaite concordance avec des solutions établies dans des conditions asymptotiques et ils sont conformes aux intuitions physiques.
© A. Albinali et al., published by IFP Energies nouvelles, 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
B : Formation volume factor (L^{3}/L^{3})
c_{t} : Total compressibility (LT^{2}/M)
D : Diffusion coefficient (L^{2}/T)
d_{w} : Fractal dimension of the walk
I_{0} : Modified Bessel functions of order zero
K : Diffusion Coefficient (L^{2}/T)
K_{0} : Modified Bessel functions of order zero
K_{1} : Modified Bessel function of order one
k_{α} : Permeability (L^{2}/T^{1−α })
k_{α,β} : Permeability (L^{ β−1}/T^{1−α })
: Matrix flux rate per unit rock volume
r_{m} : Radius of matrix sphere (L)
X_{F} : Fracture halflength (L)
Greek
Θ : Anomalous diffusion coefficient
ϑ : Bias factor for preferential diffusion direction
Subscripts
Superscripts
Abbreviations
ADMNF: Anomalous Diffusion for Matrix and Natural Fractures
CTRW: ContinuousTime Random Walk
NFR: Naturally Fractured Reservoirs
Introduction
Fluid flow in naturally fractured nanoporous reservoirs have been primarily modeled using approaches that treat the flow domain as continuum and the flow parameters as represented by their statistical averages. In these models, rock matrix and natural fractures are generally perceived as low permeability storage medium and highly conductive pathways, respectively, that exchange fluids under certain domain and flow conditions. Darcy’s law is the governing fluid flow equation with the mathematical models utilizing a normal diffusion process in the matrix and natural fracture regions — transport is considered to be random and following a Gaussiandistributed probabilitydensity function.
Normal diffusion constitutes that a particle’s mean square displacement 〈r〉 is related to time linearly; that is:(1)
where D is the diffusion coefficient. The resulting diffusivity equation used in conventional petroleum engineering fluid flow models is:(2)
Several conceptual models have been developed for Naturally Fractured Reservoirs (NFR): (i) a single medium with enhanced permeability due to natural fractures, (ii) dualporosity medium where tight rock matrix feeds into conductive natural fractures, (iii) matrix with a discrete fracture network where each fracture is individually characterized, and (iv) fractal media where properties are scaled with distance to some reference point in the domain. The current models are appropriate for conventional reservoirs, which have permeabilities ranging from tens to hundreds of milliDarcy (1 Darcy ≈ 10^{−12} m^{2}). However, the existing approaches are not adequate in describing flow in nanoporous fractured reservoirs due to the extreme petrophysical complexity and velocityfield heterogeneity. Production rates decline drastically after a short period of time when the fluid stored in the natural fractures is produced and the lowpermeability matrix cannot feed into the fractures at rates matching the fracture depletion (Nelson, 2001).
In tight, unconventional reservoirs, such as shalegas and tightoil plays, matrix and natural fracture permeabilities are downscaled significantly. The permeability of the rock matrix is in the nano to microDarcy range and the pore radii range from a few to no more than hundreds of nanometers. Natural fractures at different sizes, conductivities, and orientations constitute the discontinuities in these systems due to their large contrast to the matrix conductivity. They may be the product of sourcerock maturation, tectonic events, burial, uplifting, or the changes in the stress field. Generally, the aperture of these fractures is larger than the opening of the matrix pores and flow channels. Additionally, these fractures can be randomly distributed as data from well logs, core samples and outcrops confirm. Another level of complexity is presented in source rocks due to the presence of organic material (kerogen) with different pore structure and permeability than the inorganic matrix. (At low pressures, existence of organic material may give rise to the contribution of desorption and molecular diffusion to flow; however, in this paper, such low pressures will not be considered.) This complex environment includes multiple scales (Fig. 1), which create preferential flow paths, differences in pressure profiles, varying fluid compositions, and complex flow regimes (CamachoVelázquez et al., 2012; Ozkan, 2013).
Figure 1. Scale heterogeneity in unconventional NFR (Russian, 2013). 
From a fundamental fluid mechanics perspective, flow in nanoporous unconventional reservoirs takes place on an assembly of short and long flow paths — compared to the total medium — in addition to the heterogeneous distribution of these channels. Under these conditions, continuum mechanics, which is a major assumption of Darcy’s law, is inapplicable (Fig. 2). Diffusive flows at slower rates take place in the matrix and advection has a small contribution to flow due to the diminutive proportion of micropores. Advection, however, is the dominant transport mechanism in natural fractures. Where fractures form a network (continuum), advection in fractures dictates the global pressure distribution, which, in turn, governs the local diffusion in the matrix. Therefore, a nonlocal, hereditary transport results in the system. In light of this discussion, unconventional NFR are conceived as disordered media and flow and transport are described by anomalous diffusion and fractional calculus in this paper.
Figure 2. Knudsen number and expected flow mechanisms (modified from Roy et al., 2003). 
Several mathematical assumptions can lead to the formulation of anomalous diffusion, which range from assuming that diffusion follows a powerlaw as a function of distance, fluid particles follow a ContinuousTime Random Walk (CTRW) behavior, or that the observation and correlation scales are different (Raghavan, 2011). Anomalous diffusion relates the MeanSquare Displacement (MSD) of a particle to time by the following relation:(3)
This becomes useful in representing fluid flow in disordered, heterogeneous porous structures, specifically at varying scales. In the CTRW model, the transitioning time (or waiting time) of a particle between two jumps is described by a probability density function. Longer waiting times can be associated to an increased flux impediment in the system due to obstructions such as discontinuous fractures. As a consequence the particle MSD will grow slower than in the normal diffusion case. A laboratory study demonstrating the potential of the CTRW model to determine the degree of subdiffusion in a heterogeneous porous media is presented in Berkowitz et al. (2000). Similarly, it is possible to define stochastic particle displacement models resulting in MSD growing faster than in the normal case, hence leading to superdiffusion. The probability density function describing particle displacement is heavytailed allowing for longer particle jumps. One such model was presented by Redner (1989) for a stratified porous medium with layers of different conductivity.
Physically, sub and superdiffusion imply that the local flows are temporally and spatially convolved and cannot be described in terms of instantaneous, local gradients. The diffusion equation describing subdiffusive flow contains a time fractional derivative, while superdiffusion is described by a spacefractional derivative. As a consequence, the solutions fall in the framework of fractional calculus.
Previous attempts of using fractional calculus have utilized the concept of fractals to scale the properties spatially. For instance, O’Shaughnessy and Procaccia (1985) have modified the diffusion equation as:(4)
and defined the diffusion coefficient K in a fractal geometry by:(5)
In Equations (4) and (5), r represents the distance to some reference point in the domain, Θ being the anomalous diffusion coefficient, and D is the fractal dimension. The solution involves a steadystate conductivity, (K); that is assigned spatially and does not take into account the temporal or spatial dependence.
The aim of this paper is to discuss the construction and utility of 1D anomalous diffusion models for oil and gas production from tight, unconventional reservoirs with fractured horizontal wells. Both analytical and numerical solutions of the diffusion equation of fractional form are presented. In the analytical solution, only timefractional anomalous diffusion is considered because of the difficulties in satisfying the boundary conditions when spacefractional derivatives are involved. The numerical model takes into account both the temporal and spatial fractional derivatives as the derivatives can be approximated numerically. The focus of the analytical solution is to document and discuss the versatility of anomalous diffusion to consider various scales and connectivity of natural fractures. Numerical solution addresses two issues: handling noflow boundaries and demonstrating the combined effect of space and timefractional diffusion in pressuredistributions in the reservoir.
1 Analytical Approach
The objective in this section is to apply anomalous diffusion to cover a wider range of heterogeneities in unconventional reservoirs. The approach taken here utilizes the dualporosity idealization of fractured reservoirs but considers additional complexities: (i) natural fractures are discrete and finite in length (do not form a continuum); anomalous diffusion describes flow in the union of the matrixfracture system, (ii) there is a set of conductive and globally connected fractures where normal diffusion prevails, while the matrix includes another set of small, discontinuous fractures which gives rise to anomalous diffusion, and (iii) there are two scales of fractures; flow in globally connected fractures can be described by anomalous diffusion due to tortuous path, varying shape, nonuniform aperture, roughness, cementation, etc., and the diffusion in the matrix is also anomalous because of the existence of discontinuous microfractures.
For convenience, the analytical solution is derived in two steps. First, a transfer function is derived to describe fluid transfer from matrix to global fractures for a radialflow system (vertical well at the center of a cylindrical reservoir) shown in Figure 3. Here, the porous medium is assumed to consist of spherical matrix elements of uniform radius enveloped by natural fractures that form a connected network for flow. (The extension of the formulation to other matrix shapes is straightforward.) Pressure is uniform at the surface of the matrix elements and flux from matrix elements to the natural fractures is distributed instantaneously and uniformly over onehalf the volume of natural fractures (de Swaan, 1976; Kazemi, 1969). The transfer function is derived by considering the conditions in items (i) and (ii) above.
Figure 3. Representation of the utilized cylindrical dualporosity system (Ozkan, 2011). 
In the second step, the solution is extended to the trilinear flow model, developed by Ozkan et al. (2009) for multifractured horizontal wells in unconventional reservoirs. The trilinear model couples flow between three contiguous regions: outer reservoir, inner reservoir and hydraulic fracture (Fig. 4). The inner reservoir is considered as a dualporosity region while the outer reservoir and the hydraulic fracture are considered homogeneous (singleporosity) regions. Anomalous diffusion considerations are invoked into the trilinear flow model by using the appropriate transfer function derived in the first step for the dualporosity inner reservoir. In this paper, we only discuss the derivation of the transfer function with anomalous diffusion. Details of the solutions presented in this section can be found in Albinali (2016).
Figure 4. Flow regions in the trilinear flow model. 
1.1 Transfer Function for Anomalous Diffusion
A CTRW process (Montroll and Weiss, 1965) may be used to define the constitutive (flux) relation of timefractional anomalous diffusion in porous media. The process consists of probability functions that describe jump lengths and waiting times, and will provide the MSD in the sought power form given in Equation (3). The resulting diffusion equation will include index/indices that can be used to describe the transport behavior and the physical structure. The following equation, proposed by Metzler et al. (1994), was utilized in this paper:(6)where the diffusion exponent α is related to the diffusion index θ by:(7)
Note here that setting θ = 0, or α = 1, corresponds to normal diffusion and we recover the standard Darcy’s law. Following Raghavan (2011) and Raghavan and Chen (2013), we can write the convective flux as:(8)and the transient diffusion equation can be stated as:(9)
In general, three parameters are used to describe a fractal structure; a spectral dimension d_{s}, fractal dimension d_{f} and fractal dimension of the walk d_{w} (Mandelbort, 1983). These parameters are related by:(10)
There is no straightforward relationship between the fractal dimension and the diffusion index θ and additional information is needed to compute the fractal dimension (Chang and Yortsos, 1990; de Swaan et al., 2012). An example is given in CamachoVelázquez et al. (2008) where the CTRW concept was used following the fractional diffusion equation by Metzler et al. (1994). The authors introduced:(11)and(12)in their formulations and included data from transient and boundary dominated flow tests. Using asymptotic approximations of dimensionless equations and plotting production versus time they, were able to calculate the slopes that correspond to α and ν and compute d_{f}. Also, they showed how to relate Equations (11) and (12) to the parameters in conventional reservoir engineering models via plotting. An overview about the CTRW and diffusion on fractals can be found in Metzler and Klafter (2000). Another approach to model diffusion on fractals was presented in Barker (1988) and Chang and Yortsos (1990) by considering a fractal geometry and applying the classical conservation of mass law to obtain a diffusivity equation. The resulting models yield isotropic properties that scale with distance, similar to Equation (4) presented by O’Shaughnessy and Procaccia (1985).
In the analytical part of this work, we consider anomalous diffusion independently for both domains of the dualporosity idealization by assigning two diffusion exponents α_{m} and α_{f} for matrix and natural fractures, respectively. Deriving the equations in this manner gives flexibility in choosing normal or anomalous diffusion when extending the solution into other reservoirs. For example, in reservoirs with high matrix permeability and discontinuous fractures, the solution can be provided as normal diffusion for the matrix and anomalous diffusion for the natural fractures. Another example is for reservoirs with homogeneous and wellconnected natural fractures and heterogeneous tight matrix where one can apply normal diffusion for the natural fractures and anomalous diffusion for the matrix.
Starting with the flow in the matrix spheres, we incorporate Equation (6) into the conservation of mass equation to get the continuity equation for a slightly compressible fluid as:(13)
Note here that α_{m} is the diffusion exponent pertaining to the matrix spheres. Similarly, flow in the natural fractures is governed by:(14)where represents the flux from matrix to natural fractures and is given by:(15)
The solution of the above system of equations is obtained in the Laplace transform domain and leads to the definition of the following dualporosity transfer function similar to those defined by Barenblatt et al. (1960) and Kazemi (1969):(16)
The details of the derivation of Equation (16) is given in Appendix A. We use the transfer function for the dualporosity inner region of the trilinear model (Ozkan et al., 2009) to obtain the results discussed in the next section. An alternative procedure is given in Noetinger and Gautier (1998) and Noetinger et al. (2001). Their work also includes discretization of the developed fluid flow equation and the CTRW algorithm over the porous medium.
1.2 Verification of the Analytical Solution and Discussion of Results
Synthetic data in Table 1 are used to verify the results of the analytical solution (ADMNF) against the TriLinear Model (TLM). The models are first compared under homogeneous reservoir conditions by setting diffusion exponents of the analytical solution to unity since the trilinear model constitutes normal diffusion. Also, the dualporosity functions in both models were set to unity to represent a homogeneous reservoir. Figure 5 shows an excellent agreement between the bottomhole pressures obtained from the two models.
Figure 5. Verification for homogeneous reservoir. TLMH versus developed solution (ADMNFH). 
Data used for analytical solution verification
Next, the models are compared under dualporosity idealization. The diffusion exponents of the analytical solution are again set to unity. Several runs were conducted under different values of matrix thickness. Transmissivity (λ) and storativity (ω) ratios of de Swaan (1976) and Kazemi (1969) are used to label the trilinear model runs. Transmissivity is a ratio of the matrix flow capacity to that of the natural fractures, i.e. smaller values indicate larger matrix slabs. Storativity is a measure of the fluid storage in natural fractures to that in the matrix. These parameters are defined as:(17)and(18)where σ is a shape factor for the matrix blocks and, for the spherical geometry used in this work, σ = 12 (Kazemi et al., 1976). Solid lines correspond to the trilinear model and circles correspond to the analytical solution in Figure 6. The results show excellent agreement under the dualporosity idealization.
Figure 6. Verification under dualporosity idealization. 
After verifying the analytical solution for asymptotic cases, we use it to investigate the effects of anomalous diffusion. Figure 7 considers normal diffusion in fractures (α_{f} = 1) and anomalous diffusion in matrix for the data in Table 1. In Figure 7, the case for α_{f} = α_{m} = 1 corresponds to the conventional dualporosity model. For the values of α_{m} ≠ 1, Figure 7 displays the effect of different levels of the matrix heterogeneity. Because matrix contribution comes in at later times (after depleting the network of natural fractures), Figure 7 focuses on the bottomhole pressure at late times.
Figure 7. Results for subdiffusive flow in matrix. 
In the second case shown in Figure 8, for a fixed value of α_{m} = 1, a range of the diffusion exponent of the natural fractures, α_{f}, from 1 to 0.1 is considered. Values less than one indicate subdiffusion in the fracture network. As expected, higherpressure drops are observed as flow deviates from normal diffusion (smaller α_{f} values). This can also be interpreted that flow becomes hindered as flow deviates from normal diffusion due to increased heterogeneity of the velocity field.
Figure 8. Results for subdiffusive flow in natural fractures. 
In Figures 9 and 10, we consider some cases of mixed diffusion. Case a, shown in Figure 9, considers two combinations of normal diffusion in natural fractures with normal and subdiffusive flows in the matrix. The discrepancy between the bottomhole pressures for normal and subdiffusive matrix flows at late times (after the fracture system is depleted) indicates that the analytical model proposed in this work extends the conventional dualporosity approach to more complex matrix heterogeneities.
Figure 9. Case a for mixed diffusion process. 
Figure 10. Case b for mixed diffusion process. 
Case b shown in Figure 10 is similar to Case a in Figure 9 except for subdiffusion, instead of normal diffusion, in fractures. For practical purposes, the difference between the results for normal and anomalous diffusion is not discernible in Figure 10. This result should be expected based on the reduced conductivity of fractures in the case of anomalous diffusion; that is, even though the matrix flow capacity increases when α_{m} increases from 0.1 to 1, the limited flow capacity of fractures due to anomalous diffusion dictates how much fluid can be transferred from matrix to fractures.
To emphasize the significance of the results shown in Figures 7–10, it must be noted that the classical dualporosity approach is only capable of representing the heterogeneity caused by the contrast between the volumeaveraged properties of the matrix and fracture media. As discussed in the Introduction, when the heterogeneity reaches the limit of continuum mechanics, flow phenomena cannot be represented in terms of the average properties of the medium. Figures 7–10 indicate that the new dualporosity solution provides a convenient means of covering a wider range of heterogeneity encountered in unconventional reservoirs.
2 Numerical Approach
As noted in Introduction, the analytical difficulty in handling noflow boundaries in spacefractional anomalous diffusion makes the numerical treatment of the problem desirable. In addition, the numerical solution has the potential to extend the utility of the anomalous diffusion model to other cases of practical interest such as multiphase flow. In this section, we present an implicit finitedifference scheme to model linear anomalous diffusion for a single phase, slightly compressible fluid flowing in a finite reservoir initially at uniform pressure, with a constantrate boundary at x = 0, and a noflux boundary at x = L (Fig. 11).
Figure 11. Schematic of the linear system examined. The boundaries x = 0 and x = L correspond to the constantrate (q = constant) and noflux (q = 0) boundaries respectively. 
The fractional flux law incorporating space and time nonlocality is defined by Chen and Raghavan (2015) as:(19)where 0 < α < 1 and 0 < β < 1 are the fractional derivative exponents in time and space respectively, μ is the fluid viscosity, and k_{α,β} is the phenomenological coefficient with dimensions L^{1+β}T^{1−α}. Note that, in the asymptotic case of α = β = 1, the flux law reverts back to the wellknown Darcy’s law with k_{α,β} in dimensions of L^{2}. For α < 1 the diffusion is expected to slow down (subdiffusion) while in the case of β < 1 particle motion should accelerate (superdiffusion).
The classical mass conservation equation in 1D is given by:(20)
Combining Equations (19) and (20), and introducing the initial and boundary conditions, the following initial boundary value problem is obtained:(21) (22) (23) (24)
The spatial domain [0, L] is discretized into a blockcentered grid of I_{max} blocks and uniform block length ∆x = L/I_{max}. The grid block centers are labeled with the index x_{i} where i = 1, …, I_{max} as shown in Figure 12. The time domain [0,T] is discretized into N + 1 time steps with uniform time increment ∆t = T/N and the time steps are labeled with the index t_{n} = n∆t, n = 0, …, N (Fig. 13). The numerical approximations of the function p(x_{i}, t_{n}) are denoted by .
Figure 12. Spatial discretization of the examined system. 
Figure 13. Temporal discretization of the examined system. 
2.1 Approximation of Fractional Derivatives
The time and space integer derivatives of the mass conservation (Eq. 20) are approximated by the forward and central differences, respectively, leading to:(25)or after substituting the flux (Eq. 19):(26)where the subscripts denote the grid block interfaces and the superscripts n + 1 and n denote the simulation time steps.
TimeFractional Derivative
In Equation (26), the time fractional derivatives in the flux terms at grid block interfaces (x _{ i∓1/2}) are defined in the Caputo (1967) sense, allowing the use of integerorder initial conditions. The leftsided Caputo derivative is defined as:(27)where p − 1 < α < p, p being the smallest integer larger than α, and Γ(x) is the gamma function. Substituting Equation (27) with 0 < α < 1 into the time fractional derivatives of Equation (26) yields: (28)
Using the approach presented in Murio (2008), after rewriting Equation (28) as a summation of integrals over uniform time intervals, ∆t, and evaluating the time derivative in each integral term by the first order forward approximation, the following expression is obtained:(29)where(30)and(31)
The complete derivation of Equation (29) is presented in Appendix B and further details are presented in Holy (2016).
Equation (29) can be further rearranged to yield:(32)
Two important observations can be made from the fractional time derivative approximation in Equation (32). First, the evaluation of the derivative at t _{ n+1} requires the spacefractional pressure derivatives at all previous time steps from t _{0} to t _{ n }. Second, for α → 1, → 1 for k ≥1 and σ _{ α,∆t } → 1, and hence Equation (32) requires the evaluation of the space fractional derivative at t _{ n+1} and t _{0} only. Because the system is initially at rest, , and, thus, the equation reduces to evaluating the spatial derivative at t _{ n+1} in the asymptotic case α = 1. Hence, Equation (32) reverts to the classical implicit formulation for normal diffusion.
Figure 14 shows the weights attributed to the past spacefractional derivatives as a function of α when computing solutions for t _{ n+1} = 50 based in Equation (32). As can be seen, α values closer to one result in derivatives putting more weight on past states while α values towards zero make the derivative more local, with the asymptotic case of α = 0 reverting to the first derivative. It should also be noted that the weight coefficients naturally add up to zero for any 0 ≤ α ≤ 1 and any number of prior time steps t _{ n }. Hence, the approximation presented in Equations (29) and (32) is not a truncated series.
Figure 14. Temporal dependence of the timefractional derivative as a function of α. 
SpaceFractional Derivative
The spacefractional derivatives at the grid block interfaces are defined as weighted twosided Caputo derivatives based on the symmetric Caputo derivative presented by Klimek and Lupa (2011). The spacederivative at the interface at x _{ i+1/2} and time t _{ n+1} is:(33)where 0 < β < 1, 0 < ϑ < 1 is the weighting factor allowing to set a bias on the preferred diffusion direction on either side of the point of interest, x _{ i+1/2}. In Equation (33), and are the left and rightsided Caputo derivatives, respectively, as shown in Figure 15.
Figure 15. Left and right side contributions to the space fractional derivative computed at x _{ i+1/2}. 
The leftsided Caputo derivative (from the left boundary at x = 0 to the interface at x _{ i+1/2}) is defined as:(34)
The rightsided Caputo derivative (from the point at x _{ i+1/2} to the right boundary at x = L) is defined as in Kilbas et al. (2006):(35)
Both derivatives are discretized using the approach presented in Zhang et al. (2007), leading to the following approximations:(36)and(37)where(38)and(39)
The detailed derivation of Equations (34) and (35) is given in Appendix C and Holy (2016).
Substituting the left and rightsided Caputo approximations, Equations (34) and (35), into Equation (33), the finitedifference approximation to the spacefractional derivative at x _{ i+1/2} becomes:(40)
In a similar manner, the finitedifference approximation to the space fractional derivative at the interface at x _{ i−1/2} and time t _{ n+1} becomes:(41)
Two important observations can be made from the spacefractional derivative approximations: first, the evaluation of the twosided derivative at time t _{ n+1} requires pressures at every grid block, leading to a fully populated (I _{ max } × I _{ max }) iteration matrix; second, for for m > 1, and , the approximations in Equations (40) and (41) reduce, respectively, to:(42)and(43)which correspond to the classic central difference approximations of the first derivatives at grid block interfaces at x _{ i+1/2} and x _{ i−1/2} in case of normal diffusion. In addition, as for the time fractional derivative, the weight coefficients attributed to each grid block naturally add up to zero for any 0 < β < 1 and any number of grid blocks. Hence, the approximations presented in Equations (40) and (41) are not truncated series.
Substituting the time and spacefractional derivative approximations, Equations (29), (40) and (41), into Equation (26), an implicit finitedifference scheme for the anomalous diffusion equation is obtained at time t _{ n+1}:(44)
Treatment of Boundaries
The rate q at the constant flux boundary (x = 0) is specified explicitly in grid block 1 by:(45)
At the noflux boundary (x = L), the flux is zero and hence for grid block I _{ max } is:(46)
2.2 Verification of the Numerical Solution
The finitedifference scheme presented in the previous section is used to model the flow of a slightly compressible fluid towards a hydraulic fracture in a complex porous media at constant production rate within onefourth of the fracture drainage area. The objective is to qualitatively assess the influence of the exponents, α and β, and verify the solution. For simplicity, the fluid and rock properties are treated as constants and the synthetic data used is presented in Table 2.
Fluid and rock properties used for verification of the numerical solution
Sensitivity on TimeFractional Exponent α
Figure 16 shows the impact of the timefractional exponent α on the pressure distribution in the reservoir at different times, where β = 1, T = 200 days, N = 200, ∆t = T/N = 1 day, and the number of grid blocks, I _{ max } = 250. The case for α = 1 corresponds to normal diffusion while the cases for α < 1 correspond to subdiffusion, which is clearly shown by the smaller drainage distances at t = 200 days. In cases for α = 1 and α = 0.99, the pressure disturbance has clearly reached the outer boundary while transient flow is still dominant for α < 0.4. In addition, pressure drawdowns at the fracture face are more significant when α < 1 (Fig. 17), which is in agreement with the responses obtained by the analytical solution.
Figure 16. Pressure profiles in the studied system at times t = 10, 20, 50, 200 days as a function of α with β = 1. The case α = 1 corresponds to normal diffusion while α < 1 corresponds to subdiffusion. 
Figure 17. Pressure drop versus time at the fracture face as a function of α with β = 1. The case α = 1 corresponds to normal diffusion while α < 1 corresponds to subdiffusion. 
Sensitivity on SpaceFractional Exponent β
Figure 18 shows the impact of the space fractional exponent β on the pressure distribution in the reservoir at different times, where α = 1, T = 200 days, N = 200, ∆t = T/N = 1, the number of grid blocks, I _{ max } = 250, and the weighting factor, ϑ = 0 (for rightsided derivative contribution only). Note that in case of ϑ = 0, the iteration matrix is reduced to an upper triangular matrix. The case β = 1, corresponds to normal diffusion while the cases β < 1 correspond to superdiffusion. This is shown through the larger drainage distances as well as smaller pressure drawdowns at the fracture face (Fig. 19).
Figure 18. Pressure profiles in the studied system at times t = 10, 20, 50, 200 days as a function of β with α =1 and ϑ = 0. The case β = 1 corresponds to normal diffusion while β < 1 corresponds to superdiffusion. 
Figure 19. Pressure drop versus time at the fracture face as a function of β with α = 1 and ϑ = 0. The case β = 1 corresponds to normal diffusion while β < 1 corresponds to superdiffusion. 
From a purely qualitative point of view, the implemented finite difference scheme seems to be well suited for handling the imposed boundary conditions although neither extensive testing nor stability and convergence analysis have been performed at this stage. Another important point to address is the selection of the weighting factor ϑ. Figure 20 shows the impact of using the rightsided (ϑ = 0), symmetric (ϑ = 0.5) and leftsided (ϑ = 1) space fractional derivative at t = 10 days for different values of β, keeping all other parameters unchanged. For β = 1, the normal diffusion case is obtained and is unaffected by the weighting factor. For β < 1, although the trend of superdiffusion is observed in all three cases, the cases for ϑ = 0.5 and ϑ = 1 seem to be physically irrelevant due to the concave downward shape of the pressure profiles.
Figure 20. Pressure profiles in the studied system at t = 200 days as a function of β and ϑ. 
Conclusion
A modified flux law incorporating spatial and temporal heterogeneity of the velocity field has been coupled with the classic mass conservation equation to model anomalous diffusion in fractured nanoporous media. The model has been used to simulate the complex flow behaviors around hydraulically fractured horizontal wells in unconventional reservoirs. The resulting partial differential equation includes noninteger time and space derivatives, which are solved by means of fractional calculus. In the analytical solution only the time dependency of flux has been studied due to the difficulty of handling noflux boundary conditions when space fractional derivatives are introduced. This drawback is overcome in the numerical approach.
Sensitivities run on the timefractional exponent α are in agreement in the analytical and numerical solutions showing larger pressure drops near the flux boundary and smaller drainage distances for smaller values of α. Such cases seem to be representative of porous media dominated by tight matrix and disconnected natural fractures leading to subdiffusion. The spacefractional exponent β has the opposite effect. Decreasing β values lead to smaller pressure drops at the flux boundary and larger drainage distances at a particular time. Such cases might represent porous media dominated by a wellconnected natural fracture network leading to superdiffusion.
References
 Albinali A. (Forthcoming 2016) Analytical Modeling of Fractured NanoPorous Reservoirs, PhD Dissertation, Colorado School of Mines. [Google Scholar]
 Barenblatt G.E., Zheltov I.P., Kochina I.N. (1960) Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, Journal of Applied Mathematics and Mechanics 24, 5, 1286–1303. [CrossRef] [Google Scholar]
 Barker J. (1988) A generalized radial flow model for hydraulic tests in fractured rock, Water Resources Research 24, 10, 1796–1804. [CrossRef] [Google Scholar]
 Berkowitz B., Scher H., Silliman S.E. (2000) Anomalous transport in laboratoryscale, heterogeneous porous media, Water Resources Research 36, 1, 149–158. [CrossRef] [Google Scholar]
 CamachoVelázquez R., Fuentescruz G., VásquezCruz M. (2008) DeclineCurve Analysis of Fractured Reservoirs with Fractal Geometry, SPE Reservoir Evaluation & Engineering 11, 3, 606–619. [CrossRef] [Google Scholar]
 CamachoVelázquez R., VásquezCruz M.A., FuentesCruz G. (2012) Recent Advances in Dynamic Modeling of Naturally Fractured Reservoirs, SPE Latin American and Caribbean Petroleum Engineering Conference, Mexico City, Mexico, April 1618. [Google Scholar]
 Caputo M. (1967) Linear Models of Dissipation whose Q is almost Frequency IndependentII, Geophysical Journal 13, 5, 529–539. [NASA ADS] [CrossRef] [Google Scholar]
 Chang J., Yortsos Y.C. (1990) Pressure transient analysis of fractal reservoirs, SPE Formation Evaluation 5, 1, 31–38. [CrossRef] [Google Scholar]
 Chen C., Raghavan R. (2015) Transient flow in a linear reservoir for spacetime fractional diffusion, Journal of Petroleum Science and Engineering 128, 194–202. [CrossRef] [Google Scholar]
 de Swaan A. (1976) Analytic Solutions for Determining Naturally Fractured Reservoir Properties by Well Testing, Society of Petroleum Engineers Journal 16, 3, 117–122. [CrossRef] [Google Scholar]
 de Swaan A., CamachoVelázquez R., VásquezCruz M. (2012) Interference Tests Analysis in Fractured Formations with a Timefractional Equation, SPE Latin American and Caribbean Petroleum Engineering Conference, Mexico City, Mexico, April 1618. [Google Scholar]
 Holy R. (Forthcoming 2016) Numerical Investigation of 1D Anomalous Diffusion in Fractured Nanoporous Reservoirs, PhD Dissertation, Colorado School of Mines. [Google Scholar]
 Kazemi H. (1969) Pressure Transient Analysis of Naturally Fractured Reservoirs with Uniform Fracture Distribution, Society of Petroleum Engineers Journal 9, 4, 451–462. [CrossRef] [Google Scholar]
 Kazemi H., Merrill Jr L.S., Porterfield K.L., Zeman P.R. (1976) Numerical Simulation of WaterOil Flow in Naturally Fractured Reservoirs, Society of Petroleum Engineers Journal 16, 6, 17–326. [CrossRef] [Google Scholar]
 Kilbas A.A., Srivastava H.M., Trujillo J.J. (2006) Fractional Integrals and Fractional Derivatives, in Theory and Applications of Fractional Differential Equations, Elsevier, Amsterdam. [Google Scholar]
 Klimek M., Lupa M. (2011) On Reflection Symmetry in Fractional Mechanics, Scientific Research of the Institute of Mathematics and Computer Science 10, 1, 109–121. [Google Scholar]
 Mandelbort B.B. (1983) The Fractal Geometry of Nature, W.H. Freeman and Company, New York. [Google Scholar]
 Metzler R., Glockle W.G., Nonnenmacher T.F. (1994) Fractional Model Equation for Anomalous Diffusion, Physica A: Statistical Mechanics and its Applications 211, 1, 13–24. [CrossRef] [Google Scholar]
 Metzler R., Klafter J. (2000) The Random Walk’s Guide to Anomalous Diffusion: a Fractional Dynamics Approach, Physics reports 339, 1, 1–77. [NASA ADS] [CrossRef] [Google Scholar]
 Montroll E.W., Weiss G.H. (1965) Random Walks on Lattices II, Journal of Mathematical Physics 6, 2, 167–181. [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Murio D.A. (2008) Implicit finite difference approximation for time fractional diffusion equations, Computers and Mathematics with Applications 56, 4, 1138–1145. [CrossRef] [MathSciNet] [Google Scholar]
 Nelson R.A. (2001) Geologic Analysis of Naturally Fractured Reservoirs, 2nd edn., Gulf Professional Publishing, Woburn. [Google Scholar]
 Noetinger B., Gautier Y. (1998) Use of the FourierLaplace transform and of diagrammatical methods to interpret pumping tests in heterogeneous reservoirs, Advances in Water Resources 21, 7, 581–590. [CrossRef] [Google Scholar]
 Noetinger B., Estebenet T., Landereau P. (2001) A direct determination of the transient exchange term of fractured media using a continuous time random walk method, Transport in porous media 44, 3, 539–557. [CrossRef] [Google Scholar]
 Ozkan E., Brown M., Raghavan R., Kazemi H. (2009) Comparison of Fractured HorizontalWell Performance in Conventional and Unconventional Reservoirs, SPE Western Regional Meeting, San Jose, California, March 2426. [Google Scholar]
 Ozkan E. (2011) On NonDarcy Flow in Porous Media: Modeling Gas Slippage in Nanopores, SIAM Mathematical & Computational Issues in the Geosciences Meeting, Long Beach, California, March 2124. [Google Scholar]
 Ozkan E. (2013) A Discourse on Unconventional Reservoir Engineering – The State of the Art after a Decade, Unconventional Reservoir Engineering Project Consortium Meeting at Colorado School of Mines, Golden, Colorado, Nov. 811. [Google Scholar]
 O’Shaughnessy B., Procaccia I. (1985) Diffusion on Fractals, Physical Review A 32, 5, 3073–3083. [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Raghavan R. (2011) Fraction Derivative: Application to Transient Flow, Journal of Petroleum Science and Engineering 80, 1, 7–13. [CrossRef] [Google Scholar]
 Raghavan R., Chen C. (2013) FracturedWell Performance under Anomalous Diffusion, SPE Reservoir Evaluation & Engineering 16, 3, 237–254. [CrossRef] [Google Scholar]
 Redner S. (1989) Superdiffusive Transport Due to Random Velocity Fields, Physica D: Nonlinear Phenomena 38, 13, 287–290. [CrossRef] [MathSciNet] [Google Scholar]
 Roy S., Raju R., Chuang H.F., Cruden B.A., Meyyappan M. (2003) Modeling gas flow through microchannels and nanopores, Journal of Applied Physics 93, 8, 4870–4879. [CrossRef] [Google Scholar]
 Russian A. (2013) Anomalous Dynamics of Darcy Flow and Diffusion through Heterogeneous Media, PhD Dissertation, Universitat Politècnica de Catalunya. [Google Scholar]
 Stehfest H. (1970) Numerical Inversion of Laplace Transforms, Communications of the ACM 13, 1, 47–49. [CrossRef] [Google Scholar]
 Zhang X., Lv M., Crawford J., Young I.M. (2007) The impact of boundary on the fractional advectiondispersion equation for solute transport in soil: Defining the fractional dispersive flux with the Caputo derivatives, Advances in Water Resources 30, 5, 1205–1217. [CrossRef] [Google Scholar]
Appendix
Appendix A: Derivation of the Analytical Solution
The mass balance equation for spherical control volume is(A. 1)
The timedependent Darcy equation given as equation governs the flux ν. The right side of equation can be written as(A. 2)
For a slightly compressible fluid, the equation can be reduced to(A. 8)
Defining dimensionless variables as(A. 9) (A. 10)and(A. 11)with(A. 12)
Taking of both sides, we can write the matrix flow equation in dimensionless form as(A. 13)
Introducing(A. 14) (A. 15)and applying Laplace transformation, we get(A. 16)which have a general solution of(A. 17)
Boundary condition(A. 18)in dimensionless form is(A. 19)and in terms of w _{ mD } is(A. 20)
This yields(A. 21)and we can write(A. 22)or equally(A. 23)
Reverting the expression to (A. 24)
The second boundary condition is(A. 25)in dimensionless form and Laplace domain is(A. 26)and(A. 27)respectively. The solution for the matrix is given as(A. 28)
For the flow in the natural fractures, we repeat steps of (A.1) to (A.8) and we need to include the flux term from the matrix region to get(A. 29)
Matrix elements supply natural fractures by(A. 30)
Taking of both sides, putting the equation in dimensionless form and applying Laplace transform we get(A. 31)
We can substitute for from previous step and defining(A. 32)and(A. 33)we arrive at(A. 34)
This is a modified Bessel’s equation of order zero with general the solution(A. 35)
The outer boundary condition is(A. 36)which leads to(A. 37)since(A. 38)
At the wellbore(A. 39)which leads to(A. 40)
Solving for C_{2} and substitute in the equation, we get(A. 41)
The solution can be inverted back to the time domain using Stehfest (1970) algorithm.
Appendix B: Finite Difference Approximation to the Time Fractional Derivative
The finite difference approximation Equation (29) to the time fractional derivative Equation (28) is derived as follows:(B. 1)
Using the approach presented in Murio (2008) this can also be written as the summation of integrals over intervals ∆t:(B. 2)
Replacing the time derivative term by the first order forward approximation:(B.3)
Integration of the integral term leads to:(B. 4)
Or after shifting indices (so that the summation starts at t _{ n+1}), the Caputo approximation becomes:(B. 5)where:(B. 6)and(B. 7)
Appendix C: Finite Difference Approximation to the Space Fractional Derivative
The discretization of the space fractional potentials at grid block interfaces is based on the approach presented in Zhang et al. (2007).
The finite difference approximation Equation (36) to the leftsided space fractional derivative Equation (34) is derived as follows:(C. 1)
The integral can be rewritten as summation of integrals over grid block lengths ∆x from grid block 1 to i where the first derivative is approximated by the central difference at grid block interfaces x _{1+1/2} to x _{ i+1/2}:(C. 2) (C. 3)
After shifting indices in Equation (C.3) (so that the summation starts from node x _{ i }) the left sided Caputo approximation in space becomes:(C. 4)where(C. 5)and(C. 6)
Using the same approach, the finite difference approximation Equation (37) of the rightsided space fractional derivative Equation (35) is derived as follows:(C. 7)
The integral can be rewritten as summation of integrals over grid block lengths ∆x from grid block i to I _{ max } where the first derivative is approximated by the central difference at grid block interfaces x _{ i−1/2} to :(C. 8) (C. 9)
Or in compact formation:(C. 10)where σ _{ β,∆x } and are defined as in (C.5) and (C.6) respectively.
Cite this article as: A. Albinali, R. Holy, H. Sarak and E. Ozkan (2016). Modeling of 1D Anomalous Diffusion in Fractured Nanoporous Media, Oil Gas Sci. Technol 71, 56.
All Tables
All Figures
Figure 1. Scale heterogeneity in unconventional NFR (Russian, 2013). 

In the text 
Figure 2. Knudsen number and expected flow mechanisms (modified from Roy et al., 2003). 

In the text 
Figure 3. Representation of the utilized cylindrical dualporosity system (Ozkan, 2011). 

In the text 
Figure 4. Flow regions in the trilinear flow model. 

In the text 
Figure 5. Verification for homogeneous reservoir. TLMH versus developed solution (ADMNFH). 

In the text 
Figure 6. Verification under dualporosity idealization. 

In the text 
Figure 7. Results for subdiffusive flow in matrix. 

In the text 
Figure 8. Results for subdiffusive flow in natural fractures. 

In the text 
Figure 9. Case a for mixed diffusion process. 

In the text 
Figure 10. Case b for mixed diffusion process. 

In the text 
Figure 11. Schematic of the linear system examined. The boundaries x = 0 and x = L correspond to the constantrate (q = constant) and noflux (q = 0) boundaries respectively. 

In the text 
Figure 12. Spatial discretization of the examined system. 

In the text 
Figure 13. Temporal discretization of the examined system. 

In the text 
Figure 14. Temporal dependence of the timefractional derivative as a function of α. 

In the text 
Figure 15. Left and right side contributions to the space fractional derivative computed at x _{ i+1/2}. 

In the text 
Figure 16. Pressure profiles in the studied system at times t = 10, 20, 50, 200 days as a function of α with β = 1. The case α = 1 corresponds to normal diffusion while α < 1 corresponds to subdiffusion. 

In the text 
Figure 17. Pressure drop versus time at the fracture face as a function of α with β = 1. The case α = 1 corresponds to normal diffusion while α < 1 corresponds to subdiffusion. 

In the text 
Figure 18. Pressure profiles in the studied system at times t = 10, 20, 50, 200 days as a function of β with α =1 and ϑ = 0. The case β = 1 corresponds to normal diffusion while β < 1 corresponds to superdiffusion. 

In the text 
Figure 19. Pressure drop versus time at the fracture face as a function of β with α = 1 and ϑ = 0. The case β = 1 corresponds to normal diffusion while β < 1 corresponds to superdiffusion. 

In the text 
Figure 20. Pressure profiles in the studied system at t = 200 days as a function of β and ϑ. 

In the text 