Modeling of 1D Anomalous Diffusion in Fractured Nanoporous Media

Oil & Gas Science and Technology, Jul 2018

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 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 non-local 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.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 multi-porosité) 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 non-locale 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 PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://ogst.ifpenergiesnouvelles.fr/articles/ogst/pdf/2016/04/ogst150102.pdf

Modeling of 1D Anomalous Diffusion in Fractured Nanoporous Media

Oil & Gas Science and Technology - Rev. IFP Energies nouvelles Modeling of 1D Anomalous Diffusion in Fractured Nanoporous Media Ali Albinali 0 1 Ralf Holy 0 1 Hulya Sarak 0 1 Erdal Ozkan 0 1 0 Colorado School of Mines, Petroleum Engineering Department , 1600 Arapahoe Street, Golden, CO 80401 - USA 1 Albinali A. (Forthcoming 2016) Analytical Modeling of Fractured Nano-Porous Reservoirs, PhD Dissertation, Colorado School of Mines , USA - 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 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 non-local 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. NOMENCLATURE A B ct D df ds dw h I0 K K0 K1 k ka ka,b L p q ^qm R r rm s t u XF a b H h k l r / ϑ x i m f GREEK F D n SUPERSCRIPTS Time step Hydraulic fracture Dimensionless 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 Gaussiandistributed probability-density function. Normal diffusion constitutes that a particle’s mean square displacement hri is related to time linearly; that is: r2 ¼ Dt 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 ð1Þ ð2Þ 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 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 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. 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: r2 8 a ¼ 1 Normal diffusion > > >>> a 6¼ 1 Anomalous diffusion ta; where< > a > 1 Super diffusion > > >:> a < 1 Sub diffusion: ð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Þ a) km scale b) m scale c) cm scale d) μm scale Scale heterogeneity in unconventional NFR (Russian, 2013). and defined the diffusion coefficient K in a fractal geometry by: KðrÞ ¼ Kr H ð5Þ 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. 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 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 Representation of the utilized cylindrical dual-porosity system (Ozkan, 2011). 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). 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: where the diffusion exponent a is related to the diffusion index h by: 2 a ¼ 2 þ h 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: ð7Þ ð8Þ ð9Þ 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 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: 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: 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 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 dualporosity idealization by assigning two diffusion exponents am and af 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. ð10Þ ð11Þ ð12Þ Note here that am is the diffusion exponent pertaining to the matrix spheres. Similarly, flow in the natural fractures is governed by: where q^m represents the flux from matrix to natural fractures and is given by: ^qm ¼ 4pr2m 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): ð13Þ ð14Þ ð15Þ f s 2km ð Þ ¼ hfDkf rmD 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. 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. 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 and 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: k ¼ rX 2F kkm f x ¼ ð/ctÞf ð/ctÞf þ ð/ctÞm ð17Þ ð18Þ 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 (af ¼ 1Þ and anomalous diffusion in matrix for the data in Table 1. In Figure 7, the case for af ¼ am ¼ 1 corresponds to the conventional dual-porosity model. For the values of am 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 am ¼ 1, a range of the diffusion exponent of the natural 1.E+9 1.E+8 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 am 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. 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 1E+5 1E+6 1E+7 Time, hour fractures, af , 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 af 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. 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: ð19Þ where 0 < a < 1 and 0 < b < 1 are the fractional derivative exponents in time and space respectively, l is the fluid viscosity, and ka;b is the phenomenological coefficient with dimensions L1þbT 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 ka;b in dimensions of L2. 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: uðx; tÞ B Combining Equations (19) and (20), and introducing the initial and boundary conditions, the following initial boundary value problem is obtained: uð0; tÞ ¼ qð0; tÞ for t ¼ A 0 ðconstant rate boundaryÞ; ð20Þ ð21Þ ð22Þ uðL; tÞ ¼ ¼ 0 for t 0 ðno flux boundaryÞ ð24Þ The spatial domain [0, L] is discretized into a blockcentered grid of Imax blocks and uniform block length Dx = 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 Dt = T/N and the time steps are labeled with the index tn = nDt, n = 0, . . ., N (Fig. 13). The numerical approximations of the function p(xi, tn) are denoted by pin. 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: uðx; tÞ B 1 x hui B iþ12 hui B i 12 Pin i 12 iþ12 ! /ct Pinþ1 ¼ B t Pin ð26Þ where the subscripts i 12 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 (xi 12) are defined in the Caputo (1967) sense, allowing the use of integer-order initial conditions. The left-sided Caputo derivative is defined as: 0CDCtaf ðtÞ ¼ Cðp 1 aÞ 0 ð27Þ ¼ Cð1 1 ð1 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: 1 A where and 1 1 ra; t ¼ Cð1 þ aÞ t1 a where p 1 < a < p, p being the smallest integer larger than a, and C(x) is the gamma function. Substituting Equation (27) The complete derivation of Equation (29) is presented in Appendix B and further details are presented in Holy (2016). ð34Þ ð35Þ ð36Þ ð37Þ ð38Þ ð39Þ on the symmetric Caputo derivative presented by Klimek and Lupa (2011). The space-derivative at the interface at xi+1/2 and time tn+1 is: #ÞxCiþ12 DbL ð33Þ 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, xiþ1=2. In Equation (33), C0DCxbiþ1=2 and xCiþ1=2 DLb are the left and right-sided Caputo derivatives, respectively, as shown in Figure 15. The left-sided Caputo derivative (from the left boundary at x = 0 to the interface at xi+1/2) is defined as: Both derivatives are discretized using the approach presented in Zhang et al. (2007), leading to the following approximations: and where and i CDb 0 xiþ12 ¼ rb; x X xðmbÞ Pinþþ21 m m¼1 xCiþ12DLb ¼ m¼1 A ¼ 2 6 6 66 n ra; t66 þ X 66 k¼1 6 6 4 xðaÞ k þ xðkaþÞ1 xðnaþÞ1 3 @xb 77 @bp xi 12; tnþ1 k 77 @xb 777 ð32Þ 7 7 7 5 Two important observations can be made from the fractional time derivative approximation in Equation (32). First, the evaluation of the derivative at tnþ1 requires the spacefractional pressure derivatives at all previous time steps from t0 to tn. Second, for a ? 1, xðkaÞ ? 1 for k 1 and ra; t? 1, and hence Equation (32) requires the evaluation of the space fractional derivative at tnþ1 and t0 only. Because the system is initially at rest, @bp xi 1=2; t0 @xb ¼ 0, and, thus, the equation reduces to evaluating the spatial derivative at tnþ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ðkaÞ þ xðkaþÞ1 attributed to the past space-fractional derivatives as a function of a when computing solutions for tnþ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 tn. 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 The detailed derivation of Equations (34) and (35) is given in Appendix C and Holy (2016). Temporal dependence of the time-fractional derivative as a function of a. 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 xi+1/2 becomes: In a similar manner, the finite-difference approximation to the space fractional derivative at the interface at xi 1=2 and time tnþ1 becomes: > >>> þð1 : i # X xðmbÞ Pnþ1 iþ2 m m¼1 Imax i #Þ X xðmbÞ Pinþþm1 m¼1 Pnþ1 iþ1 m > >>> þð1 : #Þ i 1 # X xðmbÞ Pnþ1 iþ1 m m¼1 Imaxþ1 i X xðmbÞ Pnþ1 i 1þm m¼1 Pnþ1 i m Pnþ1 i 2þm 9 > > > > = > > > > ; ð41Þ Left and right side contributions to the space fractional derivative computed at xi+1/2. Two important observations can be made from the spacefractional derivative approximations: first, the evaluation of the two-sided derivative at time tnþ1 requires pressures at every grid block, leading to a fully populated ðI max I maxÞ iteration matrix; second, for b ! 1; xðmbÞ ! 0 for m > 1, and rb; x ! 1= x, the approximations in Equations (40) and (41) reduce, respectively, to: 0 t Pin > > > >>> þð1 : i # X xðmbÞ Pinþþ21 m m¼1 Imax i #Þ X xðmbÞ Pinþþm1 m¼1 9 9 1 Pinþþ11 m >>>> >> => >>>>>>>>>>>>= CCCCCCC > > > Pinþ11þm >>>; C > @bp xiþ12; tnþ1 k 1A >>>>>>>>>>>>>;> CCCCCCCCC @xb i 1 # X xðmbÞ Pinþþ11 m m¼1 Pinþm1 Imaxþ1 i #Þ X m¼1 xðmbÞ Pinþ11þm > Pinþ21þm > > > > ; @bp xi 12; tnþ1 k 1A >>>>>>>>>>>>>> CCCCCCA @xb ; C 9=> 9>>>>>>>>>>>>>>= CCCCCCCCCC > > > > > > C 1 ¼ x # Pnþ1 iþ1 # Pinþ1 Pinþ11 þ ð1 #Þ Pinþ1 Pinþ11 g ¼ Pinþ1 Pinþ11 which correspond to the classic central difference approximations of the first derivatives at grid block interfaces at xiþ1=2 and xi 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 tnþ1: ð43Þ Treatment of Boundaries The rate q at the constant flux boundary (x = 0) is specified explicitly in grid block 1 by: 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, 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. Pressure 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. 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, Imax = 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 Pressure 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. 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. times, where a ¼ 1, T ¼ 200 days, N ¼ 200, t ¼ TN ¼ 1, the number of grid blocks, Imax = 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. 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 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. 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. Barker J. (1988) A generalized radial flow model for hydraulic tests in fractured rock, Water Resources Research 24, 10, 1796-1804. Berkowitz B., Scher H., Silliman S.E. (2000) Anomalous transport in laboratory-scale, heterogeneous porous media, Water Resources Research 36, 1, 149-158. 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. 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. Caputo M. (1967) Linear Models of Dissipation whose Q is almost Frequency Independent-II, Geophysical Journal 13, 5, 529-539. Chang J., Yortsos Y.C. (1990) Pressure transient analysis of fractal reservoirs, SPE Formation Evaluation 5, 1, 31-38. 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. de Swaan A. (1976) Analytic Solutions for Determining Naturally Fractured Reservoir Properties by Well Testing, Society of Petroleum Engineers Journal 16, 3, 117-122. de Swaan A., Camacho-Velázquez R., Vásquez-Cruz M. (2012) Interference Tests Analysis in Fractured Formations with a Timefractional Equation, SPE Latin American and Caribbean Petroleum Engineering Conference, Mexico City, Mexico, April 16-18. Holy R. (Forthcoming 2016) Numerical Investigation of 1D Anomalous Diffusion in Fractured Nanoporous Reservoirs, PhD Dissertation, Colorado School of Mines. Kazemi H. (1969) Pressure Transient Analysis of Naturally Fractured Reservoirs with Uniform Fracture Distribution, Society of Petroleum Engineers Journal 9, 4, 451-462. 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. 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. 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. Mandelbort B.B. (1983) The Fractal Geometry of Nature, W.H. Freeman and Company, New York. 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. Metzler R., Klafter J. (2000) The Random Walk’s Guide to Anomalous Diffusion: a Fractional Dynamics Approach, Physics reports 339, 1, 1-77. Montroll E.W., Weiss G.H. (1965) Random Walks on Lattices. II, Journal of Mathematical Physics 6, 2, 167-181. Murio D.A. (2008) Implicit finite difference approximation for time fractional diffusion equations, Computers and Mathematics with Applications 56, 4, 1138-1145. Nelson R.A. (2001) Geologic Analysis of Naturally Fractured Reservoirs, 2nd ed., Gulf Professional Publishing, Woburn. 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. 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. 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. 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. 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. O’Shaughnessy B., Procaccia I. (1985) Diffusion on Fractals, Physical Review A 32, 5, 3073-3083. Raghavan R. (2011) Fraction Derivative: Application to Transient Flow, Journal of Petroleum Science and Engineering 80, 1, 7-13. Raghavan R., Chen C. (2013) Fractured-Well Performance under Anomalous Diffusion, SPE Reservoir Evaluation & Engineering 16, 3, 237-254. Redner S. (1989) Superdiffusive Transport Due to Random Velocity Fields, Physica D: Nonlinear Phenomena 38, 1-3, 287-290. 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. Russian A. (2013) Anomalous Dynamics of Darcy Flow and Diffusion through Heterogeneous Media, PhD Dissertation, Universitat Politècnica de Catalunya. Stehfest H. (1970) Numerical Inversion of Laplace Transforms, Communications of the ACM 13, 1, 47-49. 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. Manuscript submitted in May 2015 Manuscript accepted in May 2016 Published online in August 2016 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. 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 and pD ¼ qkBfhl p r rD ¼ X F gf t tD ¼ X 2F For a slightly compressible fluid, the equation can be reduced to Defining dimensionless variables as ðA:1Þ ðA:2Þ ðA:3Þ ðA:4Þ ðA:5Þ ðA:6Þ ðA:7Þ ðA:8Þ ðA:9Þ ðA:10Þ ðA:11Þ with Introducing and applying Laplace transformation, we get which have a general solution of ¼ X 2F gm wmDðrD; tDÞ ¼ pmDðrD; tDÞrD bm ¼ X 2F gm ~ gf X 2F am sam d2wmD dr2D bmwmD ¼ 0 Taking @@taamm 11 of both sides, we can write the matrix flow equation in dimensionless form as Boundary condition in dimensionless form is and in terms of wmD is This yields and we can write or equally wmD ¼ Aexp pffibffiffimffiffirD þ Bexp pffibffiffimffiffirD pmðr ¼ 0; tÞ ¼ finite pmDðrD ¼ 0; tDÞ ¼ finite wmDðrD ¼ 0; tDÞ ¼ 0 A ¼ B wmD ¼ Bhexp pffibffiffimffiffirD exp pffibffiffimffiffirD i wmD ¼ 2Bsinh pffibffiffimffiffirD ðA:12Þ ðA:13Þ ðA:14Þ ðA:15Þ ðA:16Þ ðA:17Þ ðA:18Þ ðA:19Þ ðA:20Þ ðA:21Þ ðA:22Þ ðA:23Þ Reverting the expression to pmD The second boundary condition is in dimensionless form and Laplace domain is and respectively. The solution for the matrix is given as 1 pmD ¼ rD wmD ¼ 2B ^qm ¼ 4pr2mhf 2 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 dpfD RD dRD 2kamX F hfDkaf rmDpffibffiffimffiffi Tanh pffibffiffimffiffirmD # 1 þ we arrive at u ¼ sf ðsÞ upfD ¼ 0 ðA:24Þ ðA:25Þ ðA:26Þ ðA:27Þ ðA:28Þ ðA:29Þ ðA:30Þ ðA:31Þ ðA:32Þ ðA:33Þ ðA:34Þ This is a modified Bessel’s equation of order zero with general the solution The outer boundary condition is which leads to since At the wellbore which leads to pfD ¼ C1I 0 pffiuffiRD þ C2K0 pffiuffiRD 1 s pffiuffiC2K1 pffiuffiRwD ¼ 1 s pfD ¼ 1 K0ðpffiuffiRDÞ spffiuffi K1ðpffiuffiRwDÞ Solving for C2 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: A; n ¼ 0; . . . ; N 1 ¼ Cð1 1 ð1 sÞ ð1 aÞds Using the approach presented in Murio (2008) this can also be written as the summation of integrals over intervals t: 1 ð1 ðB:2Þ ðA:35Þ ðA:36Þ ðA:37Þ ðA:38Þ ðA:39Þ ðA:40Þ ðA:41Þ ðB:1Þ Replacing the time derivative term by the first order forward approximation: Integration of the integral term leads to: ffi CðaÞ k¼1 ðk 1Þ t 1 Xnþ1 ¼ CðaÞ k¼1 Z k t ðk 1Þ t 1 Xnþ1 ¼ CðaÞ k¼1 1 1 Xnþ1 ¼ CðaÞ a k¼1 1 Xnþ1 ¼ Cð1 þ aÞ k¼1 1 Xnþ1 ¼ Cð1 þ aÞ k¼1 sÞ ð1 aÞds ðtnþ1 sÞ1 ð1 aÞ#k t 1 ð1 aÞ ðk 1Þ t ½ ðtnþ1 s a k t Þ ðk 1Þ t ððn þ 1Þ t k tÞa # þððn þ 1Þ t ðk 1Þ tÞa " ðn þð1n þ k1 þ 1kÞÞaa # ta ðtnþ1 ðB:3Þ @bp xi 12; tk 1 1" @xb A ðn ðn k þ 2Þa # k þ 1Þa Or after shifting indices (so that the summation starts at tnþ1), the Caputo approximation becomes: where: and 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: C0DCxbiþ12 ¼ Cð1 1 bÞ 0 Z xiþ12 @Pðn; tnþ1Þ xiþ12 @n dn; i ¼ 1; . . . ; I max 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 x1þ1=2 to xiþ1=2: x C0DCxbiþ12 ffi Cð1 1 Xi Z m x bÞ m¼1 ð m 1Þ x Pnþ1 mþ1 Pnþ1 m xiþ12 n b dn ¼ Cð1 ¼ Cð2 1 i X bÞ m¼1 1 i X bÞ m¼1 Pnþ1 mþ1 x Pnþ1 mþ1 ¼ Cð2 ¼ Cð2 1 1 i X bÞ m¼1 1 Xi bÞ xb m¼1 x Pnþ1 m Pnþ1 mþ1 Pnþ1 mþ1 x xiþ12 1 n b 1 b3m x 7 5 ð m 1Þ x Pnmþ1 h ði x m xÞ1 b þ ði x ð m 1Þ xÞ1 bi hði þ 1 mÞ1 b ði mÞ1 bi x1 b Pnþ1 hði þ 1 m mÞ1 b ði mÞ1 bi After shifting indices in Equation (C.3) (so that the summation starts from node xi) the left sided Caputo approximation in space becomes: where and i C0DCxbiþ12 ¼ rb; x X xðmbÞ Pinþþ21 m m¼1 ðm 1Þ1 b; xð1bÞ Using the same approach, the finite difference approximation Equation (37) of the right-sided space fractional derivative Equation (35) is derived as follows: xiþ12 dn; i ¼ 2; . . . ; I max ðC:1Þ ðC:2Þ ðC:3Þ ðC:4Þ ðC:5Þ ðC:6Þ ðC:7Þ The integral can be rewritten as summation of integrals over grid block lengths x from grid block i to Imax where the first derivative is approximated by the central difference at grid block interfaces xi 1=2 to xImax 1=2: n dn ðm Or in compact formation: where rb; x and xðmbÞ are defined as in (C.5) and (C.6) respectively. ðC:8Þ ðC:9Þ ðC:10Þ Pinþ11þm


This is a preview of a remote PDF: https://ogst.ifpenergiesnouvelles.fr/articles/ogst/pdf/2016/04/ogst150102.pdf

Ali Albinali, Ralf Holy, Hulya Sarak, Erdal Ozkan. Modeling of 1D Anomalous Diffusion in Fractured Nanoporous Media, Oil & Gas Science and Technology, 56, DOI: 10.2516/ogst/2016008