Variable-order derivative time fractional diffusion model for heterogeneous porous media

Author’s Accepted Manuscript Variable-order derivative time fractional diffusion model for heterogeneous porous media Abiola D. Obembe, M. Enamul Hossain, Sidqi A. Abu-Khamsin www.elsevier.com/locate/petrol

PII: DOI: Reference:

S0920-4105(17)30377-7 http://dx.doi.org/10.1016/j.petrol.2017.03.015 PETROL3903

To appear in: Journal of Petroleum Science and Engineering Received date: 16 August 2016 Revised date: 12 January 2017 Accepted date: 7 March 2017 Cite this article as: Abiola D. Obembe, M. Enamul Hossain and Sidqi A. AbuKhamsin, Variable-order derivative time fractional diffusion model for heterogeneous porous media, Journal of Petroleum Science and Engineering, http://dx.doi.org/10.1016/j.petrol.2017.03.015 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1

Variable-order derivative time fractional diffusion model for heterogeneous porous media Abiola D. Obembe a, M. Enamul Hossain a,b,*, Sidqi A. Abu-Khamsin a a b

King Fahd University of Petroleum & Minerals, Dhahran 31261, Saudi Arabia Memorial University of Newfoundland, St. John's, NL, Canada, A1B 3X5

[email protected] [email protected] [email protected] *Corresponding author.

Abstract Constant-order derivative (COD) time fractional diffusion models have been successfully employed to describe transport processes where the rate of diffusion or the diffusion behavior differs from the classical Brownian motion (i.e. anomalous diffusion). More so, the Variable-order derivative (VOD) time fractional diffusion models have been recently acknowledged as an alternative and more rigorous approach to handle time-dependent anomalous diffusion behavior observed in highly heterogeneous fractured porous media. This study exhibits two finite difference approximations based on control volume finite difference approximations to handle a VOD time fractional diffusion mathematical model describing fluid flow in a heterogeneous porous medium. Subsequently, the numerical models are validated against the exact solution of the simplified COD time fractional diffusion mathematical model with constant rock and fluid properties. The wellestablished COD time fractional diffusion mathematical model is a special case of the herein presented VOD time fractional diffusion mathematical model. Furthermore, a modified incremental balance check and a modified wellbore flow model are presented to account for the variable order fractional exponent. The results from the numerical simulations are presented to reveal the peculiar features of the VOD time fractional diffusion model. This study would herald more applications of the VOD time fractional diffusion models in numerical modelling of fluid flow in porous media of fractal geometry or heterogeneous media.

Keywords: Constant-order derivative; time fractional diffusion equation; variable-order derivative; finite difference.

2

Nomenclature Cross-sectional area of rock perpendicular to the flow of flowing fluid in x direction, ft

2

Oil formation volume factor, bbl/stb Oil compressibility, psia-1 Total compressibility of the systems, psia-1 Formation rock compressibility of the systems, psia-1 Diameter of spherical grains, ft Reservoir height, ft Control volume centroid counter Reservoir permeability, mD Pseudo-permeability, mD. Reservoir length along x direction, ft Gridblock counter for multidimensional flow Old time level New time level Pressure of the system, psia Initial pressure of the system, psia Flowing bottom hole pressure, psia Source term, stb/day Wellbore radius, ft Time, day filtration velocity in x direction, ft/day Flow dimension at any point along the x-direction, ft Greek Symbols Fractional order of differentiation, dimensionless Fractional order derivative, , dimensionless Fractional change in viscosity per unit change of pressure, psia-1 Tortuosity, dimensionless Oil density, lb/ft3 Size of gridblock in x direction, ft Reference density, lb/ft3 Porosity, fraction Initial porosity, fraction Reference porosity, fraction Proportional parameter, dimensionless Oil specific gravity fluid dynamic viscosity, cP Oil viscosity above bubble point pressure, cP Oil viscosity at bubble point pressure, cP Ratio of the pseudo-permeability of the medium with memory to fluid viscosity, mD. /cP Acronyms and Field Units American petroleum institute

3 Reservoir barrel Standard barrel Standard cubic feet Conversion factors 1 foot

0.3048 m

1 psia

6.894757 kPa

1cp

0.001 Pa.s

1bbl/day

0.1589873 std m3/day

1 lbm/ft3

16.01846 kg/m3 0.9869233E-6 m2

1 bbl/stb

1 m3/std m3

1. Introduction Subterranean hydrocarbon reservoirs consist of rock, hydrocarbon fluids, and water. The rock serves as the storage space and is known to play an important role in fluid flow in the reservoirs. The basic law describing the diffusion process of fluids through a porous medium is the well-known Darcy law. Most authors who study diffusion problems in porous media use this equation. Similarly, almost all commercial reservoir simulators employ Darcy’s equation as a constitutive equation relating the volumetric flux to the fluid potential (or pressure) gradient. Unfortunately, Darcy’s equation was formulated based on several simplified assumptions. Thus, it is necessary to understand those assumptions before its application in the fluid mass balance. Various modifications or corrections to Darcy equation have been presented in the literature accounting for different observed phenomena i.e. slip, inertia, non-isothermal flow, anomalous diffusion, and non-Darcy flow etc. (Bell and Nur, 1978; Biot, 1957; Boley and Tolins, 1962; Caputo, 1998). The literature on modelling anomalous diffusion in reservoir rock is still burgeoning as compared to its applications in other fields of science and engineering where the anomalous diffusion models have been employed with great success (Chen and Liu, 2016; Moghaddam et al., 2016) . Anomalous diffusion occurs when a particle plume spreads at a rate that cannot be predicted by classical Brownian motion models. That is, in systems where the mean square displacement (i.e. variance) grows faster or perhaps slower than the Gaussian diffusion process. A general expression between the mean square variance and time, presented in Eq. (1) was first reported in Ben-Avraham and Havlin (2000).

4

where

(1) {

The application of fractional derivatives (in time and/or space) in describing such behavior has led to the increasing growth of novel mathematical formulations using the anomalous diffusion approach via non-local constitutive temporal and spatial flux equations. In general, the main advantage of fractional calculus (i.e. fractional derivatives) is that it provides an excellent tool for the description of memory and other hereditary properties of various materials and processes (Podlubny, 1998). Numerous numerical schemes have been reported in the literature (Huang et al., 2016; Kumar et al., 2016; Sousa and Li, 2015; Zaky et al., 2016) for various onedimensional fractional diffusion models based on constant diffusion coefficient, or COD time fractional diffusion models. Scherer et al. (2008) presented an explicit/implicit finite-difference scheme applicable to the treatment of some fractional extensions of the energy balance used for predicting the temperature field in an oil stratum. The fractional derivative was interpreted in the Grünwald–Letnikov’s sense. Murio (2008) proposed an implicit finite-difference scheme for a time-fractional (Caputo derivative) diffusion equation. Liu et al. (2009) presented an implicit finite-difference scheme for a fractional sub-diffusion equation with a nonlinear source term. Sweilam et al. (2012) proposed a finite-difference scheme that could handle efficiently a class of nonlinear, fractional Riccati differential equation. Mustapha (2010) proposed and analyzed two hybrid implicit schemes to solve both fractional super-diffusion and sub-diffusion equations. The hybrid schemes consist of time-stepping Crank–Nicolson method and finite elements for the spatial discretization. Variable time steps were employed to compensate for the singular behavior of the solution. Recently, Mustapha and AlMutawa (2012) formulated a finite difference scheme using Crank–Nicolson in time combined with the second-order central finite difference for the spatial discretization. Unfortunately, practical problems do not allow for such oversimplification of constant coefficient of diffusion (linear diffusion models), or the COD time fractional diffusion models (Chechkin et al., 2005). VOD time fractional diffusion models have been recognized as suitable to characterize complex diffusion encountered in highly heterogeneous fractured or disordered porous media (Sun et al., 2012, 2009; Umarov and Steinberg, 2009). For instance, these models are comfortably used in systems where the diffusion behavior may change with time, spatial locations or system parameters (Sun et al., 2014, 2010). These models may also be employed in systems where the diffusion behavior change from anomalous diffusion to normal diffusion and vice versa (Addison et al., 1998; Sokolov et al., 2004). It is worth mentioning that the classic diffusion models with time dependent diffusion coefficients have been successfully employed to fit some experimental data, as well as in some cases, to simulate time dependent diffusion phenomena. It is our opinion that the COD time fractional diffusion equations may not capture the complexities of the flow problem in a porous media in its entirety. Thus, following the hypothesis by Sun et al. (2009), a time dependent VOD time fractional diffusion model is employed as the right approach to describe the diffusion behavior in

5 complex (heterogeneous) porous media. Other applications of the VOD time fractional diffusion models can be found in the literature (Lin et al., 2009; Samko, 1995; Samko and Ross, 1993; Zhuang et al., 2009). In general, petroleum engineers seldom have sufficient information to effectively characterize sub-surface geological formations due to the complex geology (e.g. cracks, fissures, and discontinuities) with significant variations in properties sometimes observed on several scales. To avoid these limitations, the COD and VOD anomalous diffusion model were proposed (Awotunde et al., 2016; Obembe et al., 2016; Raghavan, 2012; Sun et al., 2014). These models are advantageous because they do not require detailed information on the matrix and fracture petro-physical properties as one would require in the normal diffusion type formulations. The coefficients associated with the COD and VOD diffusion models are phenomenological and have different physical interpretations from coefficients associated with the classic normal diffusion-based models. Due to the complex nature of the variable-order anomalous diffusion models, the only technique to handle such problems involves the development of efficient numerical schemes that can cater for such difficulties (nonlinearity, variable coefficients, etc.). The main motivation of this work is twofold; first, the well-known COD non-local constitutive flux equation relating fluid flux to pressure gradient (Caputo, 1998; Raghavan, 2012, 2011) is extended or modified to a generalized VOD non-local constitutive equation. Second, two numerical schemes are presented to handle the resulting VOD time fractional diffusion mathematical model describing fluid flow in porous media. To the authors knowledge, no such work has been reported in the petroleum engineering literature or in the literature at large. The structure of this paper is organized as follows: In Section 2, a step by step derivation of the nonlinear VOD time fractional diffusion mathematical model describing fluid flow in a hydrocarbon reservoir is presented. Two finite difference approximations, a modified incremental balance check, and modified wellbore flow model are presented based on the VOD time fractional differential operator in Section 3. In section 4 the numerical scheme validation and the results of the numerical simulations are presented and discussed to demonstrate diffusion behavior through such systems. The conclusions with some useful remarks are reported in Section 5.

2. Theoretical development 2.1. Time dependent VOD time fractional diffusion model To determine the pressure distribution in a reservoir rock, the conservation of mass (mass balance), a constitutive equation and an equation of state are combined to obtain the governing equation for the single phase flow of a slightly compressible fluid. For this study, a one-dimensional Cartesian reservoir fully saturated with oil is considered as illustrated in Fig. 1.

6 If we take a block within the reservoir, and apply the mass conservation principle for the fluid entering through the left boundary at , and leaving through the right boundary i.e.

, the final form of the governing differential equation describing fluid

flow is presented in Eq. (2). Interested readers can refer to (Ertekin et al., 2001) for full details. (

)

(

)

(2)

where is the volume conversion factor, expressed as oil field units (or , in the metric unit).

5.615, for the customary

Constitutive equation To incorporate the heterogeneity at different scales in porous media of fractal geometry, a generalized time dependent non-local flux law is presented to describe fluid flow in porous media as follows ( )

⃑

(

)

(3)

( ) where is the pseudo-permeability with dimension [ L2 ], is the transmissibility conversion factor, expressed as 1.127 for the customary oil field units (or , in the metric unit), and is the oil viscosity in centipoise (cP).

Two types of VOD time fractional operator definition exist in the literature (Coimbra, 2003; Lorenzo and Hartley, 2002). The first definition, considers the derivative order has no memory related to past derivative order values (Coimbra, 2003; Sun et al., 2009), and the second definition considers the derivative order has memory related to past derivative order values (Lorenzo and Hartley, 2002). However, in the present study, the definition of the VOD operator in which the derivative order has no memory of past derivative order values is employed. Therefore, the VOD time fractional operator is expressed as follows (Sun et al., 2012): ( )

( )

[

(

∫ )]

( ) [

(4)

( )] ( )

( )

where

Substituting Eq. (3) into Eq. (2) results in:

{

( )

( )}

(

)

(5)

7 The equation of state Equation of state relates the density as a function of pressure and temperature. Since most oil field liquid systems are considered slightly compressible, the equation of state can be represented by; (6) where is the oil density at standard condition, and factor. Incorporating Eq. (6) into Eq. (5) results in: ( )

{

( )}

is the oil formation volume

( )

(7)

where For convenience, Eq. (7) can be re-written as follows: {

( )

( )

}

( )

(8)

For slightly compressible fluids, the Right Hand Side (RHS) of Eq. (8) can be expressed in conservative form as follows (Jamal et al., 2006) : ( ) where

(

)

(9)

is the total compressibility of the saturated porous medium.

Thus substituting Eq. (9) into Eq. (8) results to: {

( )

( )

The coefficient

}

(

( )

,

and

)

(10) .

Equation 10 is the final form of the VOD time fractional diffusion equation describing fluid flow in a porous medium. The proposed Eq. (10) is a non-linear time fractional PDE due to the dependence of the coefficients on both sides on the unknown pressure. Two different finite difference approximations (i.e. Explicit and Implicit based finite difference approximations) are presented to solve the Eq. (10) numerically. Please refer to nomenclature for definition and units of other terms introduced above. It is worth noting that for (i.e. approaches 1), the Caputo time fractional derivative reduces to the classic first order time derivative. (

)

(

)

(

)

(11)

8

3. Finite difference approximations The notation are introduced in this section, where is the size of the gridblock, and is the integration time-step. In the VOD time fractional diffusion, central difference is employed to approximate the second order spatial derivative. The time derivative is approximated as shown in Section 3.1 and 3.2.

3.1. Explicit finite difference discretization First, Eq. (10) is discretized using an explicit forward finite difference scheme. To ( ) do so, the variable-order Caputo type time fractional derivative operator is approximated as follows (Karatay et al., 2011; Murio, 2008; Zhuang and Liu, 2006): ( ) where,

( )

∑

(

and

)

( )

(

Thus using Eq. (12), we approximate ( as follows: {

( )

}

(

[ (

)

(

)

)]

(12)

.

) on the control volume centered at (xi, tn) by

)

(13)

Employing the central difference approximation of the first derivative in Eq. (13) leads to ( )

[{

}

{

( )

}

]

(

Now, use of central differences to approximate (

)

)

and (

)

(14)

yields

(

)

(15)

(

)

(16)

Substituting Eqs. (15) and (16) into Eq. (14) results in ( where gridblocks.

( )

) , and

(

)

(

)

(17)

is the harmonic average between two adjacent

9 The explicit difference approximation implies the base time level is equation. Thus the forward Euler approximation of Eq. (17) results in (

)

(

where the compressibility term

(

) (

, is expressed as:

in the flow )

(18)

]

(19)

)

Rearranging Eq. (18). [

(

)

(

)

Equation 19 is the discretized form of Eq. (10), where explicit finite difference method was employed. This discretized form of Eq. (10) is unique because it can handle the variable-order fractional order term. Furthermore, Eq. (19) is only applicable for , where is the number of spatial control volume. The non-linearity in Eq. (10) was handled by approximating the transmissibility at the unknown pressure with the pressure at the previous time step. This kind of linearization has been accepted to be a reasonable approximation for slightly compressible fluids due to the weak dependence of the fluid properties (i.e. the transmissibility) on fluid pressure. By extension of the linear stability analysis presented by Awotunde et al., (2016) for the COD time fractional diffusion model based on the Caputo interpretation of the fractional operator. The VOD time fractional diffusion model presented in this study is conditionally stable for: ̅

(

(20)

)

By engineering thinking or approach demonstrated by Abou-Kassem et al. (2006), Eq. (19) can be generalized for two-dimensional and full field (three-dimensional) simulation studies for single phase flow using the following shorthand notation: [∑

(

)

∑

]

(21)

where m is an index counter, is the set whose elements are the existing neighboring gridblocks in the reservoir, is the set whose elements are the reservoir boundaries that are shared by gridblock m, is the flow rate of the fictitious well representing fluid transfer between reservoir boundary l and gridblock m as a result of boundary condition.

10

3.2. Implicit finite difference discretization Furthermore, Eq. (10) is discretized using an implicit backward Euler finite difference scheme. This implies the base time level is taken at . Furthermore, the ( )

variable-order Caputo type time fractional derivative operator as follows (Karatay et al., 2011; Zhuang and Liu, 2006): (

)

( )

∑

[ (

)

will be approximated

(

)]

)

(

(22) Simplifying further results in: {[ ( where,

(

)

( )] ( )

and

)

( )

∑ (

[ (

)

)]} (23)

.

On further simplification of Eq. (23) yields: [ (

)

( )]

Thus using Eq. (24), we approximate ( as follows: {

( )

}

( )

∑

(

[ (

)

(

)]

(24)

) on the control volume centered at (xi, tn) by

)

(25)

Employing the central difference approximation of the first derivative in Eq. (25) leads to ( )

[{

}

{

( )

}

]

(

)

(26)

(

)

(27)

Substituting Eqs. (15) and (16) into Eq. (26) results in ( Similarly,

) ( )

, and

(

)

is the harmonic average between two

adjacent gridblocks. The implicit difference approximation implied in the base time level is equation. Thus the backward Euler approximation of Eq. (27) results in

in the flow

11 {[

(

)] {[

(

)}

(

(

)]

(

)}

) (28)

where ( )

∑

[

]

(29)

Rearranging Eq. (28) ( {

) [

(

Equation 30 is applicable for

)

]}

(30)

.

Equation 30 is the discretized form of Eq. (10) where implicit finite difference method is employed. This discretized equation, i.e. Eq. (30), is unique because it can handle the variable-order exponent fractional order term, and is only applicable for , where is the number of spatial control volume. Awotunde et al., (2016) established that the implicit finite difference scheme for the analogous COD time fractional diffusion equation based on the Caputo interpretation of the fractional derivative operator is unconditionally stable. Thus, by extension, the above formulated implicit finite difference scheme is unconditionally stable. For the implicit difference formulation, the transmissibility terms are handled using a simple iteration method. The nonlinearity is evaluated one iteration behind the pressure solution. A Successive Over Relaxation (SOR) iterative method was implemented to handle the non-linearity of the fractional equation. In summary, to solve the discretized equation to be presented for the implicit scheme, a guess pressure solution is assumed and then solved until the guessed pressure and the obtained pressure satisfy a convergence criterion. In this study, the convergence criterion is set as: (

|

)

( ) ( )

|

(31) Similarly, for two-dimensional or full field (three-dimensional) simulation studies for single phase flow, Eq. (30) can be re-written in short hand form as: ∑

(

)

∑

(

)

(32)

12 where the variables in Eq. (32) are as defined earlier.

3.3.Wellbore model The wellbore pressure is usually calculated using the steady-state model proposed by Peaceman (1983). However, due to the memory formalism, some modification is required to obtain wellbore pressure. Starting with Eq. (3), and considering a cylindrical wellbore located at gridblock i, the volumetric rate can be expressed as: (

)

where Gw is a geometric factor, derivative of the wellbore pressure (

is the Caputo variable-order time fractional . The geometric factor is expressed as: (34)

)

)

√(

In Eq. (34), radius,

(33)

(

) is the equivalent radius,

is the wellbore

Recall that, [ (

)

( )]

( )

∑

[ (

)

(

)]

(35)

Thus, the wellbore pressure is calculated by: [

(

)

(

)]

(36)

)

(

)

(37)

Simplifying further (

)

(

Therefore (

) (

)

(38) where Ei and Ewf are defined via Eq. (29).

3.4. Material balance check To ascertain the numerical approximations employed to discretize the VOD time fractional diffusion equation preserves material balance, an incremental material balance

13 check is carried out at the end of each time step. The final form of the modified incremental material balance check accounting for memory derived in Appendix B as Eq. (B.3) is ∑ ∑

[( ∑

{

) (

(

) ]

)

(39)

}

In the above equations, is the set whose elements are the reservoir boundaries that are shared by gridblock i, ∑ ( ) represents the effect of the fictitious well rates introduced to account for flow across the reservoir boundaries, represents the algebraic sum of all production rates through wells within the reservoir domain, represents the a fictitious well accounting for the contribution of anomalous diffusion at each gridblock in the reservoir domain, and finally the numerator represents the rate of accumulation. To ensure material balance is preserved, the material balance coefficients should be between 0.999995 to 1.000005 (Ertekin et al., 2001). For complete derivation, please refer to Appendix B.

3.5. Analytical solution of the constant order derivative time fractional diffusion equation The analytical solution to the COD time fractional one-dimensional diffusion model was derived using the Laplace transformation technique under the following assumptions: 1. Homogeneous porous medium. 2. Constant fluid and rock properties. 3. No producer or well present within the reservoir. 4. Specified flux at left boundary. 5. No flow at right boundary. 6. Constant order fractional derivative, where . Although presented in Laplace domain, the pressure profile in the reservoir is described by Eq. (40). See Appendix C for detailed derivation and definition of coefficients associated with Eq. (40). ̂

{

√

[

(√ (√

) )

]}

4. Result and discussion

(√

)

{ √

[

( √

)

} ]

(√

)

(40)

14 In this section, a validation study and two examples are presented to demonstrate the application of the presented finite difference approximations to model anomalous diffusion in a reservoir rock. In the first study (Section 4.1), a homogeneous reservoir with no producer or injector well distributed within the reservoir domain is considered. The main objective of this study is to ascertain the accuracy of the presented numerical schemes (model validation) by comparison with the analytical solution of the same problem. The variable-order fractional derivative is modelled as follows: ( )

(

)

(41)

where is a proportional parameter, is the measured time length, and is the initial fractional order. For the case of COD time fractional diffusion problem, the proportional parameter, is equal to zero. Model calibration by fitting the above numerical model with experimental data, or field data remain the best way to estimate the introduced phenomenological coefficients. Algorithms such as Levenberg–Marquardt algorithm (LMA), differential evolution, and particle swarm optimization come to mind. However, for practical applications, it may be appropriate to initially assume that the diffusion behavior does not differ greatly from Darcy flow. Thus, choosing initial estimates of and prior to model calibration is recommended. Examples 1 (i.e. Section 4.2.2) and 2 (i.e. Section 4.2.2), considers the problem of fluid flow in 1D and 2D reservoir with producer and injector wells distributed across the flow domains respectively. In these examples, the effect of the time dependent fractional order derivative term between 0.1 and 0.3 for the VOD time fractional diffusion model is compared with the analogous COD cases of fractional order derivative term equal to 0.1 and 0.3 respectively.

4.1. Validation of proposed model with its analytical solution for onedimensional flow Consider a homogeneous reservoir with a specified rate (350 stb/day) boundary condition at the left boundary. The reservoir right boundary is considered as a no-flow boundary (i.e. is closed). The reservoir dimensions are: length 5000 feet (1524 m), width 2000 feet (609. 6 m), thickness 100 feet (30.48 m). It has an average pseudopermeability of 10 mD.dayα, a porosity of 0.35 with an initial reservoir pressure of 7500 psi (5.17 ×107 Pa). The viscosity and formation volume factor are 0.8 cP (8 10-4 Pa.s) and 1.2 bbl /stb (1.2 m3/m3) respectively. The total compressibility of the system is taken to be 7 10 −6 psi −1 (1 .015 ×10 −9 Pa −1). To solve the problem with the implicit and explicit numerical schemes, the reservoir was discretized into fifty gridblocks implying that = 100 ft (30. 48 m). Likewise, the analytical solution, see Eq. (40), was inverted back to time domain using the Stehfest algorithm (Stehfest, 1970), and was compared with our numerical schemes for values of 0.05, 0.15 and 0.25 as shown in Figs. 2 to 5.

15 The plots of Figs. 2 to 4 present the pressure profile across the reservoir length after 20 days and 50 days of simulation for both numerical schemes and the analytical solution. The Fig. 5 on the other hand, illustrates the pressure history at block 10 throughout the simulation run for both numerical schemes and the analytical solution for different values of the fractional order derivative. As can be observed from the above figures, the pressure distribution predicted from both numerical schemes matches the analytical solution for all values of . Thus, it can be concluded that the above numerical schemes are accurate and can be employed for more practical problems. Due to the stability limitations of the explicit finite difference scheme, further analysis would employ the implicit finite difference formulation.

4.2. Investigation of significance of fractional exponent using proposed model The effect of the fractional exponent, , is investigated by employing the VOD time fractional diffusion model considering two examples. First in Section 4.2.1, a heterogonous one-dimensional reservoir with two producer and two injector wells located within the reservoir domain, with no-flow boundaries. Example 2 (i.e. Section 4.2.2) demonstrates the behavior of the implicit numerical scheme considering a hypothetical heterogeneous two-dimensional reservoir cross-section. Note that in both examples, to ensure single phase flow in the reservoir, the injected fluid and the produced fluid are considered to be the same which in reality is not the case.

4.2.1. Example 1: Heterogeneous 1D reservoir with distributed source/sink In this example, we consider the fluid flow problem in a fictitious reservoir, employing the reservoir geometry described in Section 4.1. The initial petro-physical properties of the reservoir rock are generated as a Gaussian field to introduce some degree of spatial heterogeneity employing MRST toolbox developed by SINTEF ICT (Lie, 2014) . More precisely, the initial porosity was generated as a field of independent normally distributed variables convolved with a Gaussian kernel. Additionally, the Carman-Kozeny (Carman, 1956, 1937) permeability porosity relationship was applied to obtain a crude approximation of the pseudo-permeability distribution via Eq. (42). (

)

(42)

Refer to Nomenclature section for definition of variables introduced above. The resulting porosity and pseudo-permeability fields are shown in Fig. 6. To incorporate the variation of the reservoir rock and fluid with pressure, the following input variables listed in Table 1 were employed based on the correlation/relationships introduced in Appendix A. Table 1 Reservoir rock and fluid parameters and boundary condition input for simulation

16 Parameter Value 0 bbl/day (0 m3 /s) 0 bbl/day (0 m3 /s) 31° (870 kg/ m3) 129 scf/stb (23 sm3/m3) 0.748 7500 psi (5.17 10-7 Pa) 3.5 10-6 psi-1 (5.08 10-10 Pa-1) 46 10-6 psi-1 (6.67 10-9 Pa-1) 200 10 days 720 days 3.5 10-6 psi-1 (5.08 10-10 Pa-1) 1 10-6 ft (1 10-6 m) 0.81 In this study, the reservoir is considered to be distributed with two producer and two injector wells as illustrated in Fig. 7. The two producer wells P-1 (located at block 4) and P-2 (located at block 25) are operated at a constant production rate of 500 stb/day (2.3 m3s-1), and 650 stb day-1 (2.99 m3s-1), respectively. Also, the injector wells I-1 (located at block 15) and I-2 (located at block 40) are operated at a specified injection rate of 300 stb/day (1.384 m3s-1) and 500 stb/day (2.3 m3s-1) respectively. The location of the producers and injectors within the reservoir are illustrated in Fig. 7. In this section, the VOD time fractional diffusion model is solved using the developed numerical scheme for three different pairs of the proportional parameter and the order of differentiation as presented in Table 2.

Table 2 Sensitivity studies for VOD simulation runs of Example 1 Run No. Run 1 Run 2 Run 3 Run 4

0.1 0.1 0.3 0.3

0.2 0 0 -0.2

For discussion, we focus on the block and wellbore pressures at the location 4 and 40 only. Figure 8 summarizes the result of the incremental material balance checks for the four simulation runs performed. The plot (Fig. 8) exhibits the incremental material balance check at each time step to be approximately equal to one, this implies that the numerical discretization scheme ensures material balance is preserved throughout the simulation run. Hence, the results from the simulation are satisfactory for interpretation. Figure 9 demonstrates the variation of block 4 and wellbore pressure with time for different fractional order variations. From Run 2 and Run 3 which describe the COD

17 cases, results show that the block and wellbore pressure decreases as the fractional order increases. In other words, as the heterogeneity of the rock matrices increases due to the presence of obstacles, cracks etc., the pressure drop in the reservoir increases. For the case of the time dependent fractional order (Runs 1 and 4) as depicted in Fig. 9, the simulation results reveal a very interesting behavior of the block and wellbore pressure. Run 1 differs from Run 4 because the diffusion process become more subdiffusive as time increases. While Run 4 implies the diffusion process becomes less subdiffusive with time. As one would expect, Run 1 initially matches the pressure trend of Run 2, however as time increases the fractional order increases and the pressure magnitude approaches that of Run 3 since the fractional order approaches the magnitude of that of Run 3. A similar trend is observed from Run 4 where the pressure profile initially follows the trend of Run 3 and as time increases, the pressure magnitudes approaches that from Run 2. The block and wellbore pressure at location 40 are depicted in Fig. 10. In addition, Fig. 10 reveals a similar trend for the pressure variation as the diffusion process becomes sub-diffusive and vice-versa. From all numerical runs, it is clear that the variable is of great significance in the diffusion of fluid in porous media. It is worth emphasizing that the fractional order derivative can also be a function of the dependent variable such as fluid temperature in studies devoted to thermal enhanced oil recovery. During thermal displacement processes, the alteration of rock and fluid properties with time can be captured more efficiently by the presented VOD time fractional approach and a corresponding heat transport model.

4.2.2. Example 2: Highly heterogeneous 2D reservoir with distributed wells The second example considers a two-dimensional heterogeneous reservoir of length 1,600 feet (487.68 m), width of 1,600 feet (487.68 m), and thickness of 50 feet (15.24 m) with distributed wells in the form of a 5-spot pattern as illustrated in Fig. 11. The 5 spot pattern is discretized into 16 ×16 gridblocks, with each grid of dimension 100 feet ×100 feet ×50 feet (30. 48 m ×30.48m×15.24m). The initial pressure of the fictitious reservoir is assumed to be at 4,500 psi (3.103×104 kPa). The reservoir has a specified flow boundary condition of 200 bbl/day (3.68 10-4 m3/s) at the western flank (i.e. ), with a no-flow boundary condition at the eastern flank (i.e. ). In the southern 4 boundary, a constant pressure of 4,200 psi (2.9×10 kPa) is maintained. The northern boundary is subjected to a specified rate of 300 bbl/day (5.52 10-4 m3/s). Other necessary fluid and rock properties data required for the simulation are listed in Table 3. Table 3 Reservoir parameters and boundary conditions for Example 2 Parameter

Value

18 300 bbl/day (5.52 10-4 m3/s) (

)

0 psi/feet (0 kPa/m) 4200psi (2.9×104 kPa) -200 bbl/day (3.68 10-4 m3/s) 129 scf/stb (23 sm3/m3) 0.748 31° (870 kg/m3) 200 16 1600 feet (487.68 m) 4500 psi (3.103×104 kPa) 3.5 10-6 psi-1 (5.076

-7

kPa-1)

46 10-6 psi-1 (6.67 10-9 Pa-1) 50 feet (15.24 m) 0.3 feet (0.09144 m) 0.5 days 20 days For analysis, the implicit finite difference scheme i.e. Eq. (32) was employed due to the greater stability it proffers as compared to the explicit finite difference scheme i.e. Eq. (21). The distributed wells within the reservoir domain locations in the , and notation, as well as a single index as presented by Ertekin et al. (2001) and Jamal et al. (2006) via Eq. (43), and the specified volumetric flow rates are presented in Table 4. (

)

(43)

Table 4 Distributed well locations with specified rates for Example 2 qsc

qsc

Eq. (42)

(stb/day)

(m3/s)

1

1

1

300

5.52

16

1

241

400

7.36

19 1

16

16

350

6.44

16

16

256

350

6.44

8

8

120

-650

-1.20

Equation (43) was derived by counting the gridblocks in the y direction (vertical direction) first (i index), followed by the x direction (h index). After the first sweep is completed, the process is repeated until the entire blocks are swept. Similar to Example 1, spatial heterogeneity is introduced via generating an initial Gaussian (normally distributed) porosity field, and a lognormally distributed initial pseudo-permeability field as shown in Figs. 12 and 13 respectively. The MRST toolbox (Lie, 2014) was also employed for generating the random and lognormal distributions employed for simulation.

The VOD time fractional diffusion model is solved to demonstrate the behavior of the gridblock pressures and well pressure considering the parameters or coefficients of the VOD model presented in Table 5.

Table 5 Sensitivity studies for VOD simulation runs of Example 2 Run No. Run 5 Run 6 Run 7 Run 8

0.1 0.1 0.3 0.3

0.2 0 0 -0.2

Recalling Eq. (32) for 2D reservoirs, the following definition applies for an interior grid block. (

),

{(

)(

)(

)(

)},

For a gridblock that falls on one or more boundary, a fictitious well (

{ } (44) ) is employed

to incorporate the fluid transfer across the reservoir boundary, and the gridblock. This implies that if any gridblock falls on two reservoir boundaries, the set , comprises of two fictitious wells that capture the fluid transfer across these boundaries and the gridblock. The general form of the fictitious well is given by (Abou Kassem et al., 2006) :

20 (

)

(45)

To address the constant flow boundary conditions specified at the west and north boundaries i.e. and in Table 3, the specified flow rates are prorated among all boundary gridblocks that share that boundary such that the fictitious well flow rates are given by: ∑

and

∑

(46)

where and refers to the transmissibility in the x and y flow directions respectively, and is the set that contains all the boundary gridblocks and share the specified boundary condition. Figure 14 demonstrates the variation of with time according to Eq. (21) for the four simulation runs. The plot demonstrates that the VOD diffusion model can describe the COD diffusion problem as well as describe a time-dependent anomalous diffusion problem. The simulation results of the modified incremental balance check for all simulation runs carried out are demonstrated in Fig. 15(a-d). Figure 15 shows the incremental balance coefficient after every time step for all simulation runs which is approximately equal to 1. The satisfactory results obtained from the incremental balance check depicted in Fig. 15 demonstrate the accuracy of the simulation results. In addition, we analyze the predicted block pressure and wellbore pressure at gridblock 1 (i.e. injector well location), and gridblock 120 (i.e. producer well location) only. Figure 16 exhibits both the block pressure variation with time (Fig. 16(a)) and the shut in pressure variation with time (Fig. 16(b)) at gridblock 1. The simulation results suggest the pressure evolution is dependent on the variation of with time. Close observation of Fig. 16 highlights the remarkable behavior of the fractional order of differentiation on the pressure trend. First, we note that at the very early time period, the effect of the magnitude is negligible due to the overlapping of the curves at time less than 2 days. The trend of the pressure evolution after days shows that the larger the magnitude of the fractional exponent, higher pressure is observed in the gridblock and the wellbore (i.e. higher pressure buildup). The simulation runs labelled Run 5 and Run 8 vary within the extremes of the COD cases tagged Run 6 and Run 7. A closer investigation of the results obtained from Run 8 depicted in Fig. 16 shows a trend reversal as compared to other cases. This reversal in slope of the pressure evolution is not unexpected because inspection of the parameters associated with Eq. (21) reveals a negative slope for Run 8. The block pressure and wellbore pressure variation

21 with time in the producer well location (gridblock 120) is exhibited in Fig. 17(a) and Fig. 17 (b) respectively. Similar trends demonstrated in Fig. 16 can be observed from the plot (e.g. the negligible effect of the fractional exponent at time less than 2 days). Second, as one would expect there is a reversal in the trend of pressure evolution in the producer location. The simulation result shows that the rate of pressure decrease in the gridblock, and pressure drawdown observed the wellbore increases as the underlying diffusion differs from the Darcy type diffusion. Again, in simulation Run 8 a reversal in the trend of the pressure decline was observed in Fig. 17 which results from the negative slope via Eq. (21). The negative slope describes a diffusion behavior which approaches Darcy type flow. Although it is not presented in this study, similar trend demonstrated in Fig. 16 was observed for the block pressure and wellbore pressure in the other injection wells.

5. Remarks and Conclusions In this study, two finite difference approximations were successfully presented to handle the VOD time fractional mathematical model describing fluid flow in in a hydrocarbon reservoir. The presented finite difference approximations comprise of interpreting the time fractional operator in the Caputo sense, and employing a second order central finite difference approximation to the spatial derivative terms. The formulated explicit finite difference scheme is not always stable, i.e. is conditionally stable. This means that for any given value of , there are choices of and for which the formulated scheme becomes unstable and thus leading to unrealistic or spurious numerical solutions. Similarly, the simulation runs establish the unconditional stability notion of the proposed fully-implicit, finite-difference, numerical scheme. In addition, the proposed model is an extension of the classic fluid flow model (i.e. Darcy based constitutive equation) used in commercial reservoir simulators when , and COD time fractional diffusion model. Furthermore, it is noted that the presented finite difference approximations can easily be modified to handle other classes of VOD time fractional diffusion equations (e.g. space dependent, concentration/temperature dependent and/or other system parameter dependent models). However, such modifications are beyond the scope of the current study. Finally, the proposed VOD time fractional diffusion equation extends the previous COD anomalous diffusion models presented in the literature (Albinali and Ozkan, 2016; Holy and Ozkan, 2016; Ozcan et al., 2014) to describe fluid flow in naturally fractured nano-porous shale gas reservoirs. The anomalous formulation is advantageous because it can capture the heterogeneity observed at different length scales in very heterogeneous and disordered media (Chen and Raghavan, 2015).

22

Acknowledgment The authors would like to acknowledge the support provided by King Abdulaziz City for Science and Technology (KACST), through the Science & Technology Unit at King Fahd University of Petroleum & Minerals (KFUPM), for funding this work through project No. 11-OIL1661-04, as part of the National Science, Technology and Innovation Plan (NSTIP).

Appendix A: Correlations for fluid and rock properties To solve Eqs. (19) or (29), we assume that the oil viscosity above bubble point pressure can be correlated by Eq. (A.1) as presented by Almehaideb (2003). [

Where

(

)]

(A.1)

is the oil viscosity at bubble point pressure obtained from below (A.2)

Likewise, other the correlations proposed by Almehaideb (2003) are employed to determine reservoir fluid properties such as; (A.3) [

(

)]

(A.4) (A.5)

(

) (

(A.6)

)

Appendix B: Incremental Material balance check In this section, the incremental material balance checks at time level n+1 incorporating the anomalous diffusion concept is derived. Employing Eq. (28), and writing Eq. (27) for each block in the system (i = 1, 2, 3... ) and then summing up all equations result in ∑ ∑

(

[∑ [( )

( ) ]

)]

∑

(∑

) (B.1)

23 The sum of all inter-block flow terms in the reservoir are expressed by the first term on the LHS of Eq. (B.1), add up to zero, while the second term on the LHS represents the M algebraic sum of all production rates through wells ∑ , those across reservoir (∑ ∑ ), and those resulting from the effect of memory ∑ .The RHS of this equation represents the sum of the accumulation terms in all blocks in the reservoir. Therefore, Eq. (B.1) becomes ∑

(∑

)

∑

[( )

( ) ]

(B.2)

Dividing Eq. (B.2) by the term on the LHS yields ∑ ∑

[(

)

(

) ]

(∑

(B.3)

)

where ∑

[

{

]}

(B.4)

Appendix C: Derivation of Analytical Solution for Constant order derivative one-dimensional fractional diffusion equation (

)

(

)

(C.1)

where Initial Condition: (

)

Boundary conditions: (

)

( )

(C.2)

( )

(C.3)

The Laplace transform of the Caputo time fractional derivative is defined below: {

( ) }

( )

( )

.

For convenience, the rock and fluid properties are considered as constant.

(C.4)

24 Taking the Laplace Transform of Eq. (C.1) with respect to time ̂

̂

[

]

(C.5)

With the coefficients defined below are all constants .

(C.6)

Therefore Eq. (C.5) can be re-written as ̂

̂

(C.7)

Eq. (C.7) is a non-homogeneous second order differential equation hence the solution would consist of the complimentary and particular solution Taking the auxiliary solution of; ̂

̂

, gives: (C.8)

√

and then, ̂

. Hence, the complimentary solution is:

(√

)

( √

)

(C.9)

where c1 and c2 are two constants to be determined from the boundary conditions. Noting that the right-hand side (RHS) of Eq. (C.7) is independent of the space variable x, so ̂ forms a particular solution of Eq. (C.7). ̂( )

(√

)

( √

)

(C.10)

Therefore, Eq. (C.10) is a general solution of Eq. (C.7). To determine c1 and c2, we take the Laplace transform of the boundary conditions represented by Eqs. (C.2) and (C.3). From Eq. (C.2) we observe ̂( )

.

(C.11)

Simplifying ̂( )

.

(C.12)

Similarly, taking the Laplace transform for the external or right boundary, i.e. Eq. (C.3) yields ̂(

)

(C.13)

25 Thus applying the boundary conditions, Eq. (C.12) can be written as: ̂( )

(√

)

(√

)

(C.14)

Therefore (C.15)

√

Now, substituting Eq. (C.13) into Eq. (C.10) ̂(

)

(√

)

(√

)

(√

)

(√

(√

)

)

(C.16)

Substituting Eq. (C.15) into Eq. (C.16) [

] (√

√

(√

)

)

(√

)

(C.17)

Simplifying Eq. (C.17) results in [

(√

)

(√

)

] (√

(√

)

)

(C.18)

Therefore (√

√

[

)

(√

(√

)

)

(C.19) ]

On further simplification we have √

[

( √

(C.20)

)

]

Consequently, √

[

√

( √

)

]

√

[

(√

)

]

(C.21)

Thus

√

[

(√ (√

) )

]

(C.22)

Finally, the pressure profile in the reservoir in Laplace space according to Eq. (C.23) can be written as:

26

̂

{

√

[

(√ (√

) )

] }

(√

)

{ √

[

( √

)

}

(√

)

(C.23)

]

Appendix D: Algorithms for numerical solutions To obtain the pressure solution at the next time level (n+1) for the explicit finite difference scheme, the calculation procedure is illustrated in the flowchart below.

For the implicit formulation, the calculation procedure is described in Figure 19.

References Addison, P.S., Qu, B., Ndumu, A.S., Pyrah, I.C., 1998. A particle tracking model for nonFickian subsurface diffusion. Math. Geol. 30, 695–716. Albinali, A., Ozkan, E., 2016. Analytical Modeling of Flow in Highly Disordered, Fractured Nano-Porous Reservoirs, in: SPE Western Regional Meeting. Society of Petroleum Engineers. Almehaideb, R.A., 2003. Improved correlations for fluid properties of UAE crude oils. Pet. Sci. Technol. 21, 1811–1831. doi:10.1081/LFT-120024563 Awotunde, A.A., Ghanam, R.A., Al-Homidan, S.S., Nasser-eddine, T., 2016. Numerical Schemes for Anomalous Diffusion of Single-Phase Fluids In Porous Media. Commun. Nonlinear Sci. Numer. Simul. Bell, M.L., Nur, A., 1978. Strength changes due to reservoir-induced pore pressure and stresses and application to Lake Oroville. J. Geophys. Res. 83, 4469. doi:10.1029/JB083iB09p04469 Ben-Avraham, D., Havlin, S., 2000. Diffusion and reactions in fractals and disordered systems. Cambridge University Press. Biot, M.A., 1957. The Elastic Coefficients of the Theory of Consolidation. J. Appl. Mech. 24, 594–601. Boley, B.A., Tolins, I.S., 1962. Transient Coupled Thermoelastic Boundary Value Problems in the Half-Space. J. Appl. Mech. 29, 637. doi:10.1115/1.3640647 Caputo, M., 1998. Diffusion of fluids in porous media with memory. Geothermics 28, 2– 19. doi:http://dx.doi.org/10.1016/S0375-6505(98)00047-9

27 Carman, P.C., 1956. Flow of gases through porous media. Academic press. Carman, P.C., 1937. Fluid flow through granular beds. Trans. Chem. Eng. 15, 150–166. Chechkin, A. V, Gorenflo, R., Sokolov, I.M., 2005. Fractional diffusion in inhomogeneous media. J. Phys. A. Math. Gen. 38, L679. Chen, C., Raghavan, R., 2015. Transient flow in a linear reservoir for space–time fractional diffusion. J. Pet. Sci. Eng. 128, 194–202. Chen, D., Liu, W., 2016. Chaotic Behavior and Its Control in a Fractional-Order Energy Demand–Supply System. J. Comput. Nonlinear Dyn. 11, 61010. Coimbra, C.F.M., 2003. Mechanics with variable-order differential operators. Ann. Phys. 12, 692–703. Ertekin, T., Abou-Kassem, J.H., King, G.R., 2001. Basic applied reservoir simulation. Society of Petroleum Engineers Richardson, TX. Holy, R.W., Ozkan, E., 2016. A Practical and Rigorous Approach for Production Data Analysis in Unconventional Wells, in: SPE Low Perm Symposium. Society of Petroleum Engineers. Huang, S., Zhang, R., Chen, D., 2016. Stability of nonlinear fractional-order time varying systems. J. Comput. Nonlinear Dyn. 11, 31007. Jamal, H., SM, F.A., M Rafiq, I., 2006. Petroleum Reservoir Simulation: A Basic Approach. Gulf Publishing Company. Karatay, İ., Bayramoğlu, Ş.R., Şahin, A., 2011. Implicit difference approximation for the time fractional heat equation with the nonlocal condition. Appl. Numer. Math. 61, 1281–1288. Kumar, D., Singh, J., Baleanu, D., 2016. Numerical Computation of a Fractional Model of Differential-Difference Equation. J. Comput. Nonlinear Dyn. 11, 61004. Lie, K.A., 2014. An introduction to reservoir simulation using MATLAB: User guide for the Matlab reservoir simulation toolbox (MRST), SINTEF ICT. Lin, R., Liu, F., Anh, V., Turner, I., 2009. Stability and convergence of a new explicit finite-difference approximation for the variable-order nonlinear fractional diffusion equation. Appl. Math. Comput. 212, 435–445. Liu, F., Yang, C., Burrage, K., 2009. Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term. J. Comput. Appl. Math. 231, 160–176. doi:http://dx.doi.org/10.1016/j.cam.2009.02.013 Lorenzo, C.F., Hartley, T.T., 2002. Variable order and distributed order fractional operators. Nonlinear Dyn. 29, 57–98. Moghaddam, B.P., Yaghoobi, S., Machado, J.A.T., 2016. An Extended PredictorCorrector Algorithm for Variable-order Fractional Delay Differential Equations. J. Comput. Nonlinear Dyn. Murio, D.A., 2008. Implicit finite difference approximation for time fractional diffusion

28 equations. Comput. Math. with Appl. 56, 1138–1145. doi:10.1016/j.camwa.2008.02.015 Mustapha, K., 2010. An implicit finite-difference time-stepping method for a subdiffusion equation, with spatial discretization by finite elements. IMA J. Numer. Anal. doi:10.1093/imanum/drp057 Mustapha, K., AlMutawa, J., 2012. A finite difference method for an anomalous subdiffusion equation, theory and applications. Numer. Algorithms 61, 525–543. doi:10.1007/s11075-012-9547-0 Obembe, A.D., Hossain, M.E., Mustapha, K., Abu-Khamsin, S., 2016. A modified memory-based mathematical model describing fluid flow in Porous Media. Comput. Math. with Appl. In press. doi:10.1016/j.camwa.2016.11.022 Ozcan, O., Sarak, H., Ozkan, E., Raghavan, R.S., 2014. A trilinear flow model for a fractured horizontal well in a fractal unconventional reservoir, in: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Peaceman, D.W., 1983. Interpretation of well-block pressures in numerical reservoir simulation with nonsquare grid blocks and anisotropic permeability. Soc. Pet. Eng. J. 23, 531–543. Podlubny, I., 1998. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications, vol. 198 of, Mathematics in Science and Engineering. Academic press. Raghavan, R., 2012. Fractional diffusion: performance of fractured wells. J. Pet. Sci. Eng. 92, 167–173. Raghavan, R., 2011. Fractional derivatives: Application to transient flow. J. Pet. Sci. Eng. 80, 7–13. Samko, S.G., 1995. Fractional integration and differentiation of variable order. Anal. Math. 21, 213–236. Samko, S.G., Ross, B., 1993. Integration and differentiation to a variable fractional order. Integr. Transform. Spec. Funct. 1, 277–300. Scherer, R., Kalla, S.L., Boyadjiev, L., Al-Saqabi, B., 2008. Numerical treatment of fractional heat equations. Appl. Numer. Math. 58, 1212–1223. doi:http://dx.doi.org/10.1016/j.apnum.2007.06.003 Sokolov, I.M., Chechkin, A. V, Klafter, J., 2004. Fractional diffusion equation for a power-law-truncated Lévy process. Phys. A Stat. Mech. its Appl. 336, 245–251. Sousa, E., Li, C., 2015. A weighted finite difference method for the fractional diffusion equation based on the Riemann–Liouville derivative. Appl. Numer. Math. 90, 22– 37. Stehfest, H., 1970. Algorithm 368: Numerical inversion of Laplace transforms [D5]. Commun. ACM 13, 47–49.

29 Sun, H., Chen, W., Chen, Y., 2009. Variable-order fractional differential operators in anomalous diffusion modeling. Phys. A Stat. Mech. its Appl. 388, 4586–4592. Sun, H., Chen, W., Li, C., Chen, Y., 2012. Finite difference schemes for variable-order time fractional diffusion equation. Int. J. Bifurc. Chaos 22, 1250085. Sun, H., Chen, W., Li, C., Chen, Y., 2010. Fractional differential models for anomalous diffusion. Phys. A Stat. Mech. its Appl. 389, 2719–2724. Sun, H., Zhang, Y., Chen, W., Reeves, D.M., 2014. Use of a variable-index fractionalderivative model to capture transient dispersion in heterogeneous media. J. Contam. Hydrol. 157, 47–58. Sweilam, N.H., Khader, M.M., Mahdy, A.M.S., 2012. Numerical studies for solving fractional Riccati differential equation. Appl. Appl. Math. 7, 595–608. Umarov, S., Steinberg, S., 2009. Variable order differential equations and diffusion processes with changing modes. arXiv Prepr. arXiv0903.2524. Zaky, M.A., Ezz-Eldien, S.S., Doha, E.H., Tenreiro Machado, J.A., Bhrawy, A.H., 2016. An Efficient Operational Matrix Technique for Multidimensional Variable-Order Time Fractional Diffusion Equations. J. Comput. Nonlinear Dyn. 11, 61002. Zhuang, P., Liu, F., 2006. Implicit difference approximation for the time fractional diffusion equation. J. Appl. Math. Comput. 22, 87–99. Zhuang, P., Liu, F., Anh, V., Turner, I., 2009. Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term. SIAM J. Numer. Anal. 47, 1760–1781.

30 Fig. 1. Block i as a reservoir volume element in 1D flow Fig. 2. Numerical scheme validation for fractional order of 0.05 after (a) 20 days, and (b) 50 days Fig. 3. Numerical scheme validation for fractional order of 0.15 after (a) 20 days, and (b) 50 days Fig. 4. Numerical scheme validation for fractional order of 0.25 after (a) 20 days, and (b) 50 days Fig. 5. Numerical scheme validation for gridblock 10 for equals (a) 0.05, (b) 0.15, and (c) 0.25 Fig. 6. Example 1- Gaussian generated rock (a) initial porosity field (b) pseudo-permeability field Fig. 7. Example 1- 1D-reservoir schematic illustrating location of wells. Fig. 8. Example 1-Incremental material balance checks for (a) Run 1, (b) Run 2, (c) Run 3, and (d) Run 4 Fig. 9. Example 1-Example 1-(a) Block pressure, (b) Wellbore pressure at block location 4 Fig. 10. Example 1-(a) Block pressure, (b) Wellbore pressure at block location 40 Fig. 11. Example 2- Schematic of 2D- reservoir domain Fig. 12. Example 2- A 16 16 initial porosity field generated as a Gaussian field Fig. 13. Example 2- Generated 16×16 initial pseudo-permeability field with a lognormal permeability distribution Fig. 14. Example 2-Variation of fractional order of differentiation with time for simulation runs Fig. 15. Example 2-Incremental material balance checks for (a) Run 5, (b) Run 6, (c) Run 7, and (d) Run 8 Fig. 16. Example 2- (a) Block pressure variation, (b) Wellbore pressure variation, with time at gridblock 1 Fig. 17. Example 2- (a) Block pressure variation, (b) Wellbore pressure variation, with time at gridblock 120 Fig. 18. Schematic Strategy for the Explicit VOD fractional diffusion equation Fig. 19. Schematic Strategy for the Implicit VOD fractional diffusion equation

Highlights   

Presents explicit and implicit formulations for the variable-order anomalous diffusion in heterogeneous porous media for single phase flow. Presents examples to study the effect of the complex variable-order diffusion behavior on the block and wellbore pressure. Presents a modified incremental material balance equation due to the effect of rock memory to ensure the numerical scheme preserves material balance.

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49