Modeling of 1D Anomalous Diffusion in Fractured Nanoporous Media

— Fractured nanoporous reservoirs include multi-scale and discontinuous fractures coupled with a complex nanoporous matrix. Such systems cannot be described by the conventional dual-porosity (or multi-porosity) idealizations due to the presence of different ﬂ ow 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 ﬂ ow 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 non-local and hereditary nature of ﬂ ow, 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.


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 regionstransport is considered to be random and following a Gaussiandistributed probability-density function.
Normal diffusion constitutes that a particle's mean square displacement r h i is related to time linearly; that is: where D is the diffusion coefficient.The resulting diffusivity equation used in conventional petroleum engineering fluid flow models is: Several conceptual models have been developed for Naturally Fractured Reservoirs (NFR): (i) a single medium with enhanced permeability due to natural fractures, (ii) dual-porosity 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 milli-Darcy (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 velocity-field heterogeneity.Production rates decline drastically after a short period of time when the fluid stored in the natural fractures is produced and the low-permeability matrix cannot feed into the fractures at rates matching the fracture depletion (Nelson, 2001).
In tight, unconventional reservoirs, such as shale-gas and tight-oil plays, matrix and natural fracture permeabilities are downscaled significantly.The permeability of the rock matrix is in the nano-to micro-Darcy 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 source-rock 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 (Camacho-Velázquez et al., 2012;Ozkan, 2013).
From a fundamental fluid mechanics perspective, flow in nanoporous unconventional reservoirs takes place on an assembly of short and long flow pathscompared to the total mediumin 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 non-local, 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.
Several mathematical assumptions can lead to the formulation of anomalous diffusion, which range from assuming that diffusion follows a power-law as a function of distance, fluid particles follow a Continuous-Time Random Walk (CTRW) behavior, or that the observation and correlation scales are different (Raghavan, 2011).Anomalous diffusion relates the Mean-Square Displacement (MSD) of a particle to time by the following relation: 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 super-diffusion.The probability density function describing particle displacement is heavy-tailed 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 super-diffusion 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 super-diffusion is described by a space-fractional 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: and defined the diffusion coefficient K in a fractal geometry by: In Equations ( 4) and ( 5), r represents the distance to some reference point in the domain, H being the anomalous diffusion coefficient, and D is the fractal dimension.The solution involves a steady-state 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 time-fractional anomalous diffusion is considered because of the difficulties in satisfying the boundary conditions when space-fractional 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 no-flow boundaries and demonstrating the combined effect of space-and time-fractional diffusion in pressure-distributions in the reservoir.

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 matrix-fracture 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 radial-flow 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 Scale heterogeneity in unconventional NFR (Russian, 2013).
Figure 2 Knudsen number and expected flow mechanisms (modified from Roy et al., 2003).
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 one-half 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.
In the second step, the solution is extended to the tri-linear flow model, developed by Ozkan et al. (2009) for multifractured horizontal wells in unconventional reservoirs.The tri-linear model couples flow between three contiguous regions: outer reservoir, inner reservoir and hydraulic fracture (Fig. 4).The inner reservoir is considered as a dual-porosity region while the outer reservoir and the hydraulic fracture are considered homogeneous (singleporosity) regions.Anomalous diffusion considerations are invoked into the tri-linear 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).

Transfer Function for Anomalous Diffusion
A CTRW process (Montroll and Weiss, 1965) may be used to define the constitutive (flux) relation of time-fractional 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: where the diffusion exponent a is related to the diffusion index h by: Note here that setting h = 0, or a = 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: and the transient diffusion equation can be stated as: Figure 4 Flow regions in the tri-linear flow model.
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: There is no straightforward relationship between the fractal dimension and the diffusion index h and additional information is needed to compute the fractal dimension (Chang and Yortsos, 1990;de Swaan et al., 2012).An example is given in Camacho-Velázquez et al. (2008) where the CTRW concept was used following the fractional diffusion equation by Metzler et al. (1994).The authors introduced: and 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 a and m 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 a m and a 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 well-connected 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: Note here that a m is the diffusion exponent pertaining to the matrix spheres.Similarly, flow in the natural fractures is governed by: where qm represents the flux from matrix to natural fractures and is given by: The solution of the above system of equations is obtained in the Laplace transform domain and leads to the definition of the following dual-porosity transfer function similar to those defined by Barenblatt et al. (1960) and Kazemi (1969): 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 tri-linear 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.

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 Tri-Linear Model (TLM).The models are first compared under homogeneous reservoir conditions by setting diffusion exponents of the analytical solution to unity since the tri-linear model constitutes normal diffusion.Also, the dual-porosity 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.
Next, the models are compared under dual-porosity idealization.The diffusion exponents of the analytical solution are again set to unity.Several runs were conducted under different values of matrix thickness.Transmissivity (k) and storativity (x) ratios of de Swaan (1976) and Kazemi (1969) are used to label the tri-linear 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: where r is a shape factor for the matrix blocks and, for the spherical geometry used in this work, r ¼ 12 (Kazemi et al., 1976).Solid lines correspond to the tri-linear model and circles correspond to the analytical solution in Figure 6.
The results show excellent agreement under the 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 (a f ¼ 1Þ and anomalous diffusion in matrix for the data in Table 1.In Figure 7, the case for a f ¼ a m ¼ 1 corresponds to the conventional dual-porosity model.For the values of a m 6 ¼ 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.
In the second case shown in Figure 8, for a fixed value of a m ¼ 1, a range of the diffusion exponent of the natural Verification for homogeneous reservoir.TLM-H versus developed solution (ADMNF-H).As expected, higher-pressure drops are observed as flow deviates from normal diffusion (smaller a 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.
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 dual-porosity approach to more complex matrix heterogeneities.
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 a 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 volume-averaged 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 dual-porosity solution provides a convenient means of covering a wider range of heterogeneity encountered in unconventional reservoirs.

NUMERICAL APPROACH
As noted in Introduction, the analytical difficulty in handling no-flow boundaries in space-fractional 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 finite-difference scheme to model linear anomalous diffusion for a single phase, slightly compressible fluid flowing in a finite reservoir initially at Results for subdiffusive flow in matrix.
Results for subdiffusive flow in natural fractures.Verification under dual-porosity idealization.
uniform pressure, with a constant-rate boundary at x = 0, and a no-flux boundary at x = L (Fig. 11).
The fractional flux law incorporating space and time non-locality is defined by Chen and Raghavan (2015) as: where 0 < a < 1 and 0 < b < 1 are the fractional derivative exponents in time and space respectively, l is the fluid viscosity, and k a;b is the phenomenological coefficient with dimensions L 1þb T 1Àa .Note that, in the asymptotic case of a ¼ b ¼ 1, the flux law reverts back to the well-known Darcy's law with k a;b in dimensions of L 2 .For a < 1 the diffusion is expected to slow down (subdiffusion) while in the case of b < 1 particle motion should accelerate (superdiffusion).
The classical mass conservation equation in 1D is given by: Combining Equations ( 19) and ( 20), and introducing the initial and boundary conditions, the following initial boundary value problem is obtained: Schematic of the linear system examined.The boundaries x = 0 and x = L correspond to the constant-rate (q = constant) and noflux (q = 0) boundaries respectively.
Pressure Drop, Pa Time, hour Case b for mixed diffusion process.Case a for mixed diffusion process.
The spatial domain [0, L] is discretized into a blockcentered grid of I max blocks and uniform block length Dx = 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 Dt = T/N and the time steps are labeled with the index t n = nDt, n = 0, . .., N (Fig. 13).The numerical approximations of the function p(x i , t n ) are denoted by p n i .

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: or after substituting the flux (Eq.19): where the subscripts iÇ 1 2 denote the grid block interfaces and the superscripts n þ 1 and n denote the simulation time steps.

Time-Fractional 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 integer-order initial conditions.The left-sided Caputo derivative is defined as: where p À 1 < a < p, p being the smallest integer larger than a, and C(x) is the gamma function.Substituting Equation ( 27) with 0 < a < 1 into the time fractional derivatives of Equation ( 26) yields: Using the approach presented in Murio (2008), after rewriting Equation (28) as a summation of integrals over uniform time intervals, Dt, and evaluating the time derivative in each integral term by the first order forward approximation, the following expression is obtained: where and 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: 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 a ? 1, x a ð Þ k ? 1 for k !1 and r a;Á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, @ b p x iÇ1=2 ; t 0 À Á @x b ¼ 0, and, thus, the equation reduces to evaluating the spatial derivative at t nþ1 in the asymptotic case a = 1.Hence, Equation (32) reverts to the classical implicit formulation for normal diffusion.
Figure 14 shows the weights Àx attributed to the past space-fractional derivatives as a function of a when computing solutions for t nþ1 ¼ 50 based in Equation (32).As can be seen, a values closer to one result in derivatives putting more weight on past states while a values towards zero make the derivative more local, with the asymptotic case of a = 0 reverting to the first derivative.It should also be noted that the weight coefficients naturally add up to zero for any 0 a 1 and any number of prior time steps t n .Hence, the approximation presented in Equations ( 29) and ( 32) is not a truncated series.

Space-Fractional Derivative
The space-fractional derivatives at the grid block interfaces are defined as weighted two-sided Caputo derivatives based on the symmetric Caputo derivative presented by Klimek and Lupa (2011).The space-derivative at the interface at x i+1/2 and time t n+1 is: where 0 < b < 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 The left-sided Caputo derivative (from the left boundary at x = 0 to the interface at x i+1/2 ) is defined as: The right-sided 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): Both derivatives are discretized using the approach presented in Zhang et al. (2007), leading to the following approximations: and where and The detailed derivation of Equations ( 34) and ( 35) is given in Appendix C and Holy (2016).Substituting the left and right-sided Caputo approximations, Equations ( 34) and ( 35), into Equation ( 33), the finite-difference approximation to the space-fractional derivative at x i+1/2 becomes: In a similar manner, the finite-difference approximation to the space fractional derivative at the interface at x iÀ1=2 and time t nþ1 becomes: Two important observations can be made from the spacefractional derivative approximations: first, the evaluation of the two-sided derivative at time t nþ1 requires pressures at every grid block, leading to a fully populated and r b; Áx !1=Áx, the approximations in Equations ( 40) and ( 41) reduce, respectively, to: and 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 < b < 1 and any number of grid blocks.Hence, the approximations presented in Equations ( 40) and ( 41) are not truncated series.Substituting the time-and space-fractional derivative approximations, Equations ( 29), ( 40) and ( 41), into Equation ( 26), an implicit finite-difference scheme for the anomalous diffusion equation is obtained at time t nþ1 : Figure 15 Left and right side contributions to the space fractional derivative computed at x i+1/2 .

Treatment of Boundaries
The rate q at the constant flux boundary (x = 0) is specified explicitly in grid block 1 by: At the no-flux boundary (x ¼ L), the flux is zero and hence for grid block I max is:

Verification of the Numerical Solution
The finite-difference 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 one-fourth of the fracture drainage area.The objective is to qualitatively assess the influence of the exponents, a and b, 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. Pressure profiles in the studied system at times t = 10, 20, 50, 200 days as a function of a with b = 1.The case a = 1 corresponds to normal diffusion while a < 1 corresponds to subdiffusion.
Sensitivity on Time-Fractional Exponent a Figure 16 shows the impact of the time-fractional exponent a on the pressure distribution in the reservoir at different times, where b ¼ 1, T ¼ 200 days, N ¼ 200, Át ¼ TN ¼ 1 day, and the number of grid blocks, I max = 250.The case for a ¼ 1 corresponds to normal diffusion while the cases for a < 1 correspond to subdiffusion, which is clearly shown by the smaller drainage distances at t ¼ 200 days.In cases for a ¼ 1 and a ¼ 0:99, the pressure disturbance has clearly reached the outer boundary while transient flow is still dominant for a < 0:4.In addition, pressure drawdowns at the fracture face are more significant when a < 1 (Fig. 17), which is in agreement with the responses obtained by the analytical solution.
Sensitivity on Space-Fractional Exponent b Figure 18 shows the impact of the space fractional exponent b on the pressure distribution in the reservoir at different  times, where a ¼ 1, the number of grid blocks, I max = 250, and the weighting factor, # ¼ 0 (for right-sided derivative contribution only).Note that in case of # ¼ 0, the iteration matrix is reduced to an upper triangular matrix.The case b ¼ 1, corresponds to normal diffusion while the cases b < 1 correspond to superdiffusion.This is shown through the larger drainage distances as well as smaller pressure drawdowns at the fracture face (Fig. 19).
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 right-sided (# ¼ 0), symmetric (# ¼ 0:5) and left-sided (# ¼ 1) space fractional derivative at t = 10 days for different values of b, keeping all other parameters unchanged.For b ¼ 1, the normal diffusion case is obtained and is unaffected by the weighting factor.For b < 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.to the difficulty of handling no-flux boundary conditions when space fractional derivatives are introduced.This drawback is overcome in the numerical approach.
Sensitivities run on the time-fractional exponent a 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 a.Such cases seem to be representative of porous media dominated by tight matrix and disconnected natural fractures leading to subdiffusion.The space-fractional exponent b has the opposite effect.Decreasing b 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 well-connected natural fracture network leading to superdiffusion.

APPENDIX Appendix A: Derivation of the Analytical Solution
The mass balance equation for spherical control volume is The time-dependent Darcy equation given as equation governs the flux m.The right side of equation can be written as We have We define For a slightly compressible fluid, the equation can be reduced to Defining dimensionless variables as Taking @ am À1 @t am À1 of both sides, we can write the matrix flow equation in dimensionless form as Reverting the expression to " p mD The second boundary condition is in dimensionless form and Laplace domain is respectively.The solution for the matrix is given as 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 Matrix elements supply natural fractures by Taking @ a f À1 @t a f À1 of both sides, putting the equation in dimensionless form and applying Laplace transform we get We can substitute for d" p mD dr D r D ¼r mD from previous step and defining This is a modified Bessel's equation of order zero with general the solution The outer boundary condition is At the wellbore Solving for C 2 and substitute in the equation, we get 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: Using the approach presented in Murio (2008) this can also be written as the summation of integrals over intervals Át: Replacing the time derivative term by the first order forward approximation: 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: @ 1Àa @t 1Àa @ b p x iÇ 1 2 ; t n @x b 0 @ 1 A ¼ r a;Át X nþ1 k¼1 x a ð Þ k @ b p x iÇ 1 2 ; t nþ2Àk @x b À @ b p x iÇ 1 2 ; t nþ1Àk @x b 0 @ 1 A ðB:5Þ where: 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 left-sided space fractional derivative Equation ( 34) is derived as follows: x iþ 1 2 À n Àb dn; i ¼ 1; . . .; I max À 1 ð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 : where and Using the same approach, the finite difference approximation Equation (37) of the right-sided space fractional derivative Equation ( 35) is derived as follows: x iþ 1 2 @P n; t nþ1 ð Þ @n n À x iþ 1 2 Àb dn; i ¼ 2; . . .; I max ð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 x ImaxÀ1=2 : Or in compact formation: Figure 1 Figure 7 Figure 6

Figure 12 Spatial
Figure 12Spatial discretization of the examined system.
and right-sided Caputo derivatives, respectively, as shown in Figure15.

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

Figure 14 Temporal
Figure 14Temporal dependence of the time-fractional derivative as a function of a.

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

Figure 18 Pressure
Figure 18Pressure profiles in the studied system at times t = 10, 20, 50, 200 days as a function of b with a =1 and # = 0.The case b = 1 corresponds to normal diffusion while b < 1 corresponds to superdiffusion.
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 non-integer 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

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

Figure 20 Pressure
Figure 20Pressure profiles in the studied system at t = 200 days as a function of b and #.
After shifting indices in Equation (C.3) (so that the summation starts from node x i ) the left sided Caputo approximation in space becomes:

m
are defined as in (C.5) and (C.6) respectively.

TABLE 1
f , from 1 to 0.1 is considered.Values less than one indicate subdiffusion in the fracture network.