Fractional radial-cylindrical diffusivity model for levels of heterogeneity in petroleum reservoirs

The generalized Taylor's Formula is used to derive a fractional radial-cylindrical diffusivity model via a fractional conservation of mass in radial geometry in a petroleum reservoir. The result is a space-fractional generalization of the diffusivity model with an arbitrary order that explains the degree of heterogeneity of the medium (a continuous spectrum of zero to one where zero heterogeneity equates to homogeneity). An implicit unconditionally stable numerical difference scheme of the linear form of the derived model is used to illustrate deviations of the model from the homogeneous medium. The variation of the order produces additional localized pressure drops during flow which is congruous with fluid inhibition effect of heterogeneity..


Introduction
The theory of fractional calculus since its formulation in 1695 has in the last few decades been widely applied to a myriad of real life phenomena (Li et al. 2012). Fractional calculus is a branch of applied mathematics that generalizes classical differentiation and integration to calculus of arbitrary (non-integer) order. This concept practically extends the domain of the order of traditional models and thus allows the obtainment of more accurate results in modeling systems of complex structure. The advantage of using non-integer models lies in their ability to capture memory and hereditary processes of different substances (Podlubny 1998), a characteristic not exhibited by the classical integer order models.
Modeling with fractional calculus has proved to be both convincing and novel with extensive applications in several areas including physics (Barkai et al. 2000), systems biology (Yuste and Lindenberg 2001), biochemistry (Yuste et al. 2006), hydrology (Benson et al. 2000) and finance (Raberto et al. 2002). Some specific contributions of fractional calculus in modeling include the fractional modeling of heat transfer in heterogeneous media with the assumption that the heat is dispersed in the air around the beam (Sierociuk et al. 2013), anomalous diffusion studies in physics where particles spread at a rate inconsistent with the classical Brownian motion model (Pagnini et al. 2013) and applications in porous media transport among others.
The modeling of fluid flow in porous media within the petroleum reservoir features two basic equations: (1) the material balance or continuity equation and (2) the equation of motion (Darcy's Law). The combination of these equations usually with additional relationships results in partial differential equations of integer order. The continuity equation derived via the continuum approach relates the mass accumulation to the divergence (Bear and Cheng 2010). The divergence in this case is obtainable in the form of a first-order Taylor series representation with the assumption that the mass flux is linear. This assumption is however not appropriate for accounting for flux changes in large heterogeneous porous media (Wheatcraft and Meerschaert 2008). The problem occurs principally because of the variation of properties spatially (Dawe 2004;Wheatcraft and Meerschaert 2008). Hydrocarbon reservoirs are typical examples of porous media which are heterogeneous in nature, yet the effects of heterogeneity are often not well represented and this becomes only obvious late during water production before the time predicted (Boyle et al. 2000) thus sparking several studies into the heterogeneity of hydrocarbon reservoirs.
A study by Wheatcraft and Meerschaert (2008) showed that by using the generalized Taylor series for approximating the mass flux, the limitations of the classical Taylor series are puzzled out thus resulting in the fractional conservation of mass model in rectilinear geometry. In this paper, we construct a radial-cylindrical diffusivity model for near wellbore single phase analysis which explains well the effect of heterogeneity in the hydrocarbon reservoir.
The paper is organized as follows. In ''Fractional calculus'' section is a presentation of some important concepts in fractional calculus. In ''Mathematical modeling of single phase flow'' section, the proposed model is presented with emphasis on the limitations of the existing model. We conclude with a presentation of the numerical stability, convergence and simulation of the proposed model and discussion of its novelty.

Fractional calculus
There are several definitions formulated for fractional derivatives and integrals of which the most familiar are the Riemann-Liouville, Caputo and Grunwald-Letnikov definitions which customarily emerge from fractional order integrals, derivatives and differences respectively (Monje et al. 2010). The Riemann-Liouville and Caputo definitions and their differences landmarks their fields of application. Engineers for example prefer the Caputo derivatives primarily due to its uniformity of results with the integer order differential equations, whereas Mathematicians and Physicists usually apply the Riemann-Liouville (Podlubny 1998). In this paper, we use the Caputo definition for its advantage of ensuring that initial conditions of fractional differential equations are cognate to those of integer order differential equations (Caputo 2012) and uniformity with the use of Generalized Taylor's series.
Definition 1 The Gamma function CðxÞ is defined as and is convergent for x [ 0.
Definition 2 The ath Caputo fractional derivative of the function f(t) is defined as where n À 1\a\n such that n 2 Z þ .

Definition 3
The Fractional Taylor's Series is the generalization of the classical Taylor Series. For a as the order of differentiation, suppose D ka a f ðxÞ 2 Cða; b for k ¼ 0; 1; . . .; m þ 1 where 0\a 1, the Fractional Taylor Series is defined as where a n x 8x 2 ða; b (Odibat and Shawagfeh 2007).

Mathematical modeling of single phase flow Introduction
The governing equations for a single phase flow evolves from a combination of fundamental concepts i.e., the conversation of mass, conservation of momentum and an equation of state for the porous medium when considered as a control volume. The shape of the control volume with which the system is modeled depends on the coordinate system which should conform as closely as possible to the flow geometry. The form of geometries in petroleum reservoir modeling is either rectilinear or curvilinear (Turgay et al. 2001), the latter, features the cylindrical coordinate exclusively used in single-well simulation problems (Ahmed 2006). Based on the cylindrical geometry displayed in Fig. 1, the control volume with sides Dr, Dh and Dz and with initial point ðr; h; zÞ, the mass flux through the r-(1265) face is where 1 c = volume conservation factor, q = fluid density and u r = volumetric velocity in the r-direction. The mass outflow at the r þ Dr end of the control volume is obtained by summing the mass flux through the r-face and the change that take place in the r-direction multiplied by the length Dr over which the change acts. This results in The relation in (5) is the Taylor Series expansion about the point r of the first order of which the exact expansion is given by Equation (6) is the classical Taylor series expansion about the point r and the assumption of using (5) as a form of approximation for (6) is that changes in the mass flux are linear or piecewise linear if considered as nonlinear.
The expression representing the difference between the mass entering the volume at r and the mass leaving at r þ Dr is given by where the bulk volume of the control element, Similar to the derivation for the mass exchange in the radial component within the element, the mass exchange in the h-direction is given by Also the z-directional component is equally given by By equating the sum of the divergence in r; h and z coordinates to the mass accumulation term, i.e., V b oð/qÞ=ot, and dividing through the equation by 1 c V b , the result is the mass-balance equation of the form where / = porosity. Equation (10) is the three dimensional radial-cylindrical mass-balance equation without external sink/source. This is because the well is assumed to be at the center of the drainage area and the external conditions are incorporated in the form of boundary conditions (Turgay et al. 2001).

Fractional radial-cylindrical diffusivity model
In order to remove the restriction that the flux has to be linear and the control volume has to be infinitesimal (Wheatcraft and Meerschaert 2008), we replace the classical Taylor series approximation for the divergence in each coordinate r; h and z with the fractional first-order Taylor series approximation of (3). The component of the flux through the r þ Dr section of the radial-cylindrical control volume is given by This implies that the net flux in the r-direction is given by The net mass flux in the h and z directions are given by (13) and (14) respectively. For an anisotropic medium, the fractional orders are not the same in all directions. Therefore, we have associate a h -order and a z -order to the h-direction and z-direction.
The radial fractional mass-balance model formulated by equating the net mass flux through the control volume to the mass accumulation term is obtained as Equation (15) represents the three dimensional radialcylindrical fractional conservation of mass or fractional continuity equation for a heterogeneous but anisotropic medium without external source/sink terms. It is important to note here that when the fractional orders are individually equal to one, i.e., a r ¼ a h ¼ a z ¼ 1, we obtain the classical three dimensional conservation of mass equation for radial geometry.
The second phase in the construction of the fluid flow model for the petroleum reservoir involves the introduction of conservation of momentum (Darcy's law) in radialcylindrical coordinates. The relations for the superficialvelocity components in the r; h, and z directions assuming that the potential gradient is equivalent to the pressure gradient are given by where b c -transmissibility conversion factor, k i -permeability in the direction of i, p-fluid pressure and l-fluid viscosity.
By substituting the expressions for the fluid components in (16) into (15) and setting the fluid density, q ¼ q sc =B where q sc is the fluid density at standard conditions and B is the formation volume factor, we obtain We now turn our attention to the mass accumulation term on the right-hand side of (17). We assume here that there is only vertical compressibility in the control volume, i.e. the control volume changes over time only in the vertical dimension, Dz. We can then express the mass accumulation term as o ot The first term in (18) which has time derivative of the change in the vertical dimension is expressed as where b s is the coefficient of compressibility for the porous medium same as obtained by Wheatcraft and Meerschaert (2008) for linear geometry. Also we consider the second term in (18) The third term in (18) involving the time derivative of the porosity is similarly simplified using the relation / ¼ / o ½1 þ c R ðp À p o Þ where / o is the porosity at standard conditions and c R is the rock compressibility. Therefore The mass accumulation term can be simplified by substituting (19), (20) and (21) into (18) Since the compressibility, c is small (10 À5 À 10 À6 psi À1 ) (Turgay et al. 2001) and c R is equally small, the mass accumulation can be expressed as where c t is the total compressibility of the fluid and the formation, i.e. c t ¼ c þ c R / o =/.
Substituting (23) into (17) we have 1 r Equation (24) is the fractional radial-cylindrical diffusivity model for a heterogeneous medium with allowance for vertical compressibility of the control volume with each coordinate having an independent order. The classical radial-cylindrical diffusivity model with vertical compressibility is obtained if the fractional orders are equal to 1, i.e. a r ¼ a h ¼ a z ¼ 1 as The fractional radial-cylindrical diffusivity model in one dimension (r direction) for a slightly compressible liquid in a heterogeneous medium for a r ¼ a is given by For a constant permeability reservoir, k r = constant, l is constant for slightly compressible liquids. We replace where For small compressibility c, ½1 þ cðp À p 0 Þ % 1 and the nonlinear term cðo a p=or a Þðop=orÞ is negligible and thus the linear fractional diffusivity model in one-dimensional radial geometry is 1 r Evaluating the Caputo derivatives and simplifying (28), we obtain 1 r a Dr aÀ1 Cða þ 1ÞCð2 À aÞ

Stability
Theorem 1 The implicit numerical system defined by the linear difference equation in (35) has a unique and unconditionally stable solution.
We prove the uniqueness of the solution via the application of the Gerschgorin theorem, every eigenvalue of the matrix A ¼ fa ij g associated with (35) defined by has a magnitude greater than unity.
Let k be an eigenvalue of the matrix A. Based on the Grunwald weights, then for 1\c 2, the results indicate that g j ! 0 for j ! 2. Also for any m, If z ¼ À1 is substituted into (38), it can be seen that P 1 m¼0 g n ¼ 0 with g 1 ¼ À1 being the only negative term in the sequence of the Grunwald weights. This means that Àg 1 ¼ c is the largest estimate and Àg 1 ! P k m¼0;m6 ¼1 g m for any k ¼ 0; 1; 2; . . .. Based on the Gerschgorin theorem, all eigenvalues of A ¼ fa ij g are contained in n disks, The a ii 's are given by a ii ¼ W 3;i À g 1 W 2;i þ 1 ¼ W 3;2 þ cW 2;1 þ 1 for i ¼ 1 a ii ¼ Àg 1 W 2;i þ 1 ¼ cW 2;1 þ 1 for i 6 ¼ 1 and the radii s i 's are computed by Since W ij ! 0 and from both cases of a ii , it implies that every eigenvalue k of the matrix has a real part great than 1. Therefore, the spectral radius of the inverse matrix qðAÞ 1. This means that the numerical scheme has unique solution.
We use the fractional Von Neumann stability procedure to show the unconditional stability of (35) by letting p k i ¼ n n e Iqih where I ¼ ffiffiffiffiffiffi ffi À1 p , q real. By substituting this into (35), we obtain W 1;i e Iqðiþ1Þh À e IqðiÀ1Þh n nþ1 À W 2;i X iþ1 j¼0 g j n nþ1 e IqðiÀjþ1Þh À W 3;i n nþ1 e Iqh À W 4;i n nþ1 e Iq0h À n nþ1 e Iqih ¼ Àn n e Iqih By dividing through (39) by e Iqih and further simplification, we obtain By substituting Euler's formula e Ih ¼ cos h AE I sin h, the amplification factor l max ¼ jn nþ1 j=jn n j is given by Using the properties of the grundwald weights, that is, P 1 k¼0 g j ¼ 0 and P m k¼0 g j \0, it means that the denominator in (41) is simplified to W 2 3;i cos 2 qði À 1Þh þ 2W 3;i ðW 4;i cos qði À 1Þh þ 1Þ þðW 4;i cos qði À 1Þh þ 1Þ 2 . Therefore l max 1 and hence the scheme is unconditionally stable.

Numerical example
In this section, a numerical example of the proposed model in one dimension is presented. We consider a reservoir with properties as given in Table 1.
In addition we use a logarithmically spaced radial coordinate with wellbore radius r w = 0.25 ft and a drainage radius r e = 2500 ft. A constant initial pressure pðr; t ¼ 0Þ ¼ 4000 psi and a Neumann problem using the specified flow rate is considered. The results indicate as shown theoretical that there is a perfect correlation between the two models.
The decreasing variation of a registering a lowering response of p wf is analogous to that of a positive skin in petroleum reservoirs due to artificial drilling occurring few feet from the wellbore as shown in Fig 4. A further consideration of the deviations of the proposed model is presented in Fig 5 with relatively strong and sharp response near the wellbore compared to the weak response further away from the wellbore. Table 2 shows a summary of statistics for the bottom hole response to a.
The results indicate that although pressure drops are experienced in all cases where a\1, these pressure drops are registered more materially as a ! 0. This phenomenon is absent in classical models analogous to the proposed model when a ¼ 1 in which case the model exhibits uniformity. For increasing levels of heterogeneity, the reservoir experiences more variability thus resulting in the pressure drop experienced near the wellbore.

Conclusions
The proposed fractional space radial-cylindrical diffusivity model can be used to capture scale-dependent effects of real reservoirs which is important in applications. In addition, the proposed model exhibits fractional divergence, accounts for resistance to flow and can be used to account for the level of heterogeneity associated with each unique reservoir. The results are similar to the skin factor due to well damage. Heterogeneity has been attributed to the inhibition of fluids (Sharp et al. 2003) and thus results l max ¼ 1 ÀW 2;i P iþ1 j¼0 g j cos qðj À 1Þh þ W 3;i cos qði À 1Þh þ W 4;i cos qih þ 1 ð41Þ n nþ1 ¼ n n ÀW 1;i e Iqh À e ÀIqh ð ÞÀW 2;i P iþ1 j¼0 g j e Iqð1ÀjÞh þ W 3;i e Iqð1ÀiÞh À W 4;i e Iiqh þ 1 ð40Þ