Dossier: Characterisation and Modeling of Low Permeability Media and Nanoporous Materials
Open Access
Issue
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 71, Number 4, Juillet–Août 2016
Dossier: Characterisation and Modeling of Low Permeability Media and Nanoporous Materials
Article Number 56
Number of page(s) 25
DOI https://doi.org/10.2516/ogst/2016008
Published online 09 August 2016

© A. Albinali et al., published by IFP Energies nouvelles, 2016

Licence Creative CommonsThis 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

A : Area (L2)

B : Formation volume factor (L3/L3)

ct : Total compressibility (LT2/M)

D : Diffusion coefficient (L2/T)

df : Fractal dimension

ds : Spectral dimension

dw : Fractal dimension of the walk

h : Height (L)

I0 : Modified Bessel functions of order zero

K : Diffusion Coefficient (L2/T)

K0 : Modified Bessel functions of order zero

K1 : Modified Bessel function of order one

k : Permeability (L2)

kα : Permeability (L2/T1−α )

kα,β : Permeability (L β−1/T1−α )

L : Length (L)

p : Pressure (M/L/T2)

q : Flow rate (L3/T)

: Matrix flux rate per unit rock volume

R : Radial distance (L)

r : Spherical distance (L)

rm : Radius of matrix sphere (L)

s : Laplace parameter

t : Time (T)

u : Velocity (L/T)

XF : Fracture half-length (L)

Greek

α : Time fractional exponent

β : Space fractional exponent

Θ : Anomalous diffusion coefficient

θ : Diffusion index

λ : Mobility ratio

μ : Fluid viscosity (M/L/T)

σ : Shape factor

ϕ : Porosity (L3/L3)

ϑ : Bias factor for preferential diffusion direction

ω : Storativity ratio

Subscripts

i : Coordinate

m : Matrix

f : Natural fractures

F : Hydraulic fracture

D : Dimensionless

Superscripts

n : Time step

Abbreviations

ADMNF: Anomalous Diffusion for Matrix and Natural Fractures

CTRW: Continuous-Time Random Walk

MSD: Mean-Square Displacement

NFR: Naturally Fractured Reservoirs

TLM: Tri-Linear Model

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 Gaussian-distributed probability-density 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) 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 m2). 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 nano-meters. 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).

thumbnail 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 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.

thumbnail 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 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:(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 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:(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 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.

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 dual-porosity 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 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.

thumbnail Figure 3.

Representation of the utilized cylindrical dual-porosity system (Ozkan, 2011).

In the second step, the solution is extended to the tri-linear flow model, developed by Ozkan et al. (2009) for multi-fractured 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 (single-porosity) 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 dual-porosity 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).

thumbnail Figure 4.

Flow regions in the tri-linear 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 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:(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 ds, fractal dimension df and fractal dimension of the walk dw (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 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:(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 df. 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 dual-porosity 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 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:(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 dual-porosity 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 dual-porosity 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.

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 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.

thumbnail Figure 5.

Verification for homogeneous reservoir. TLM-H versus developed solution (ADMNF-H).

Table 1.

Data used for analytical solution verification

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 (λ) and storativity (ω) 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:(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 tri-linear model and circles correspond to the analytical solution in Figure 6. The results show excellent agreement under the dual-porosity idealization.

thumbnail Figure 6.

Verification under dual-porosity 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 dual-porosity 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.

thumbnail 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, higher-pressure 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.

thumbnail 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 dual-porosity approach to more complex matrix heterogeneities.

thumbnail Figure 9.

Case a for mixed diffusion process.

thumbnail 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 710, it must be noted that the classical dual-porosity 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 710 indicate that the new dual-porosity 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 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 uniform pressure, with a constant-rate boundary at x = 0, and a no-flux boundary at x = L (Fig. 11).

thumbnail Figure 11.

Schematic of the linear system examined. The boundaries x = 0 and x = L correspond to the constant-rate (q = constant) and no-flux (q = 0) boundaries respectively.

The fractional flux law incorporating space and time non-locality 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 L1+βT1−α. Note that, in the asymptotic case of α = β = 1, the flux law reverts back to the well-known Darcy’s law with kα,β in dimensions of L2. 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 block-centered grid of Imax blocks and uniform block length ∆x = L/Imax. The grid block centers are labeled with the index xi where i = 1, …, Imax 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 tn = n∆t, n = 0, …, N (Fig. 13). The numerical approximations of the function p(xi, tn) are denoted by .

thumbnail Figure 12.

Spatial discretization of the examined system.

thumbnail 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.

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:(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 space-fractional 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 space-fractional 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.

thumbnail Figure 14.

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

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:(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 right-sided Caputo derivatives, respectively, as shown in Figure 15.

thumbnail Figure 15.

Left and right side contributions to the space fractional derivative computed at x i+1/2.

The left-sided Caputo derivative (from the left boundary at x = 0 to the interface at x i+1/2) is defined as:(34)

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):(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 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:(40)

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:(41)

Two important observations can be made from the space-fractional 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 (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 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:(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 no-flux boundary (x = L), the flux is zero and hence for grid block I max is:(46)

2.2 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, α 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.

Table 2.

Fluid and rock properties used for verification of the numerical solution

Sensitivity on Time-Fractional Exponent α

Figure 16 shows the impact of the time-fractional 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.

thumbnail 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.

thumbnail 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 Space-Fractional 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 right-sided 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).

thumbnail 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.

thumbnail 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 right-sided (ϑ = 0), symmetric (ϑ = 0.5) and left-sided (ϑ = 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.

thumbnail 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 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 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 α 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 space-fractional 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 well-connected natural fracture network leading to superdiffusion.

References

  • Albinali A. (Forthcoming 2016) Analytical Modeling of Fractured Nano-Porous 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. [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 laboratory-scale, heterogeneous porous media, Water Resources Research 36, 1, 149–158. [CrossRef] [Google Scholar]
  • Camacho-Velázquez R., Fuentes-cruz G., Vásquez-Cruz M. (2008) Decline-Curve Analysis of Fractured Reservoirs with Fractal Geometry, SPE Reservoir Evaluation & Engineering 11, 3, 606–619. [Google Scholar]
  • Camacho-Velázquez R., Vásquez-Cruz M.A., Fuentes-Cruz G. (2012) Recent Advances in Dynamic Modeling of Naturally Fractured Reservoirs, SPE Latin American and Caribbean Petroleum Engineering Conference, Mexico City, Mexico, April 16-18. [Google Scholar]
  • Caputo M. (1967) Linear Models of Dissipation whose Q is almost Frequency Independent-II, Geophysical Journal 13, 5, 529–539. [Google Scholar]
  • Chang J., Yortsos Y.C. (1990) Pressure transient analysis of fractal reservoirs, SPE Formation Evaluation 5, 1, 31–38. [Google Scholar]
  • Chen C., Raghavan R. (2015) Transient flow in a linear reservoir for space-time 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. [Google Scholar]
  • de Swaan A., Camacho-Velázquez R., Vásquez-Cruz M. (2012) Interference Tests Analysis in Fractured Formations with a Time-fractional Equation, SPE Latin American and Caribbean Petroleum Engineering Conference, Mexico City, Mexico, April 16-18. [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. [Google Scholar]
  • Kazemi H., Merrill Jr L.S., Porterfield K.L., Zeman P.R. (1976) Numerical Simulation of Water-Oil Flow in Naturally Fractured Reservoirs, Society of Petroleum Engineers Journal 16, 6, 17–326. [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. [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. [Google Scholar]
  • Montroll E.W., Weiss G.H. (1965) Random Walks on Lattices II, Journal of Mathematical Physics 6, 2, 167–181. [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 Fourier-Laplace 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. [Google Scholar]
  • Ozkan E., Brown M., Raghavan R., Kazemi H. (2009) Comparison of Fractured Horizontal-Well Performance in Conventional and Unconventional Reservoirs, SPE Western Regional Meeting, San Jose, California, March 24-26. [Google Scholar]
  • Ozkan E. (2011) On Non-Darcy Flow in Porous Media: Modeling Gas Slippage in Nano-pores, SIAM Mathematical & Computational Issues in the Geosciences Meeting, Long Beach, California, March 21-24. [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. 8-11. [Google Scholar]
  • O’Shaughnessy B., Procaccia I. (1985) Diffusion on Fractals, Physical Review A 32, 5, 3073–3083. [Google Scholar]
  • Raghavan R. (2011) Fraction Derivative: Application to Transient Flow, Journal of Petroleum Science and Engineering 80, 1, 7–13. [Google Scholar]
  • Raghavan R., Chen C. (2013) Fractured-Well Performance under Anomalous Diffusion, SPE Reservoir Evaluation & Engineering 16, 3, 237–254. [Google Scholar]
  • Redner S. (1989) Superdiffusive Transport Due to Random Velocity Fields, Physica D: Nonlinear Phenomena 38, 1-3, 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. [Google Scholar]
  • Zhang X., Lv M., Crawford J., Young I.M. (2007) The impact of boundary on the fractional advection-dispersion 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 time-dependent Darcy equation given as equation governs the flux ν. The right side of equation can be written as(A. 2)

We have(A. 3)and(A. 4)

We define(A. 5)

Now(A. 6)

Then(A. 7)

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 C2 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 left-sided 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 right-sided 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

Table 1.

Data used for analytical solution verification

Table 2.

Fluid and rock properties used for verification of the numerical solution

All Figures

thumbnail Figure 1.

Scale heterogeneity in unconventional NFR (Russian, 2013).

In the text
thumbnail Figure 2.

Knudsen number and expected flow mechanisms (modified from Roy et al., 2003).

In the text
thumbnail Figure 3.

Representation of the utilized cylindrical dual-porosity system (Ozkan, 2011).

In the text
thumbnail Figure 4.

Flow regions in the tri-linear flow model.

In the text
thumbnail Figure 5.

Verification for homogeneous reservoir. TLM-H versus developed solution (ADMNF-H).

In the text
thumbnail Figure 6.

Verification under dual-porosity idealization.

In the text
thumbnail Figure 7.

Results for subdiffusive flow in matrix.

In the text
thumbnail Figure 8.

Results for subdiffusive flow in natural fractures.

In the text
thumbnail Figure 9.

Case a for mixed diffusion process.

In the text
thumbnail Figure 10.

Case b for mixed diffusion process.

In the text
thumbnail Figure 11.

Schematic of the linear system examined. The boundaries x = 0 and x = L correspond to the constant-rate (q = constant) and no-flux (q = 0) boundaries respectively.

In the text
thumbnail Figure 12.

Spatial discretization of the examined system.

In the text
thumbnail Figure 13.

Temporal discretization of the examined system.

In the text
thumbnail Figure 14.

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

In the text
thumbnail Figure 15.

Left and right side contributions to the space fractional derivative computed at x i+1/2.

In the text
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Figure 20.

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

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.