Numerical Fractional-Calculus Model for Two-Phase Flow in Fractured Media

Numerical simulation of two-phase flow in fractured porous media is an important topic in the subsurface flow, environmental problems, and petroleum reservoir engineering.The conventionalmodel does notworkwell inmany cases since it lacks thememory property of fracture media. In this paper, we develop a new numerical formulation with fractional time derivative for two-phase flow in fractured porous media. In the proposed formulation, the different fractional time derivatives are applied to fracture and matrix regions since they have different memory properties. We further develop a two-level time discrete method, which uses a large time step for the pressure and a small time step size for the saturation.The pressure equation is solved implicitly in each large time step, while the saturation is updated by an explicit fractional time scheme in each time substep. Finally, the numerical tests are carried out to demonstrate the effectiveness of the proposed numerical model.


Introduction
Numerical simulations for multiphase flow in fractured porous media are very important in the subsurface flow, environmental problems, and petroleum reservoir engineering.Compared to usual heterogenous porous media, the fractured media have two spaces with two distinct scales: the fracture and the matrix.The fractures possess the higher permeability than the matrix, but their volume is very small when compared to the matrix.Several conceptual models [1][2][3][4][5][6][7][8][9][10][11] have been developed to simulate the multiphase flow in porous media, for example, the single-porosity model, the dual-porosity/dual-permeability model, and the discretefracture model.The review for these models can be found in [4].
These models can deal well with the two-phase flow with constant physical parameters.However, the fluid flow may perturb the porous formation by causing particle migration resulting in pore clogging or chemically reacting with the medium to enlarge the pores or diminish the size of the pores [12].As a matter of fact, the porosity and permeability often depend on the fluid pressure and saturation.For incompressible two-phase flow in petroleum reservoir engineering, the porosity varies small with respect to the pressure.However, the injection fluids can change the porosity to a great extent.In the displacement process of oil by water, the injected water will make the media wet, and then the porosity will be enlarged, along with the change of permeability.As a result, the variability of the physical parameters should be considered in modeling realistic two-phase flow.For flow in porous media, Caputo [13] has studied the behavior of fluxes in porous media using a memory formalism, in which the ordinary time derivative is replaced by a fractional derivative.Garra [14] has studied a fractional time derivative generalization of a previous Natale-Salusti model about nonlinear temperature and pressure waves, propagating in fluid-saturated porous rocks.
In this paper, we will study a fractional time derivative generalization of a classical two-phase flow model.The ordinary time derivative is replaced by the Caputo fractional derivative with variable lower limit of integral.Because the Caputo fractional derivative is composed of the convolution of a power law kernel and the ordinary derivative of the function, it is a useful instrument to describe a power law frequency variability of the physical coefficients.Since the porosity and permeability have a remarkable contrast Advances in Mathematical Physics between matrixes and fractures, the proposed formulation uses the different fractional time derivatives for the flow in fracture and matrix regions to represent their different memory properties.
Many numerical methods have been developed in the literature, for example, [15][16][17].In order to simulate the fractional two-phase flow, we propose a two-level time discrete method based on the physical property that the pressure varies less rapidly than the saturation [18]; that is, a large time step is used for the pressure, along with a small time step size for the saturation.The local fractional differential equations have been widely studied in [19].In our method, the memory of saturations is restricted within each large time step since the pressures are determined by saturations without historical memory.The proposed method is also viewed as the fractional generalization of the classical IMplicit Pressure Explicit Saturation (IMPES) method [20,21], which is a popular time-stepping approach employed in multiphase flow simulation.Like IMPES, we split the coupled system into one pressure equation and one saturation equation based on the property of multiphysics processes of two-phase flow and treat the saturation and capillary pressure in the pressure equation explicitly to eliminate its nonlinearity.Each large time step is further divided into a few substeps, in which the saturation is updated by an explicit fractional time scheme.The cell-centered finite difference method [22] is employed for spatial discretization.Finally, numerical results are given to demonstrate the validity of the proposed numerical model.

Fractional Model of Two-Phase Incompressible Flow
2.1.Fractional Model.The Caputo fractional derivative is given by the convolution of a power law kernel and the ordinary derivative of the function.So it is a useful instrument to consider a power law frequency variability of the coefficients by a simple convolution.We now introduce the definition of the Caputo fractional derivative [17] as where () is a function and 0 <  < 1.Here, we use the modified Caputo fractional derivative; that is, the lower limit of integral in (1) is taken to be a function of  instead of the constant as usual; that is,  = ().The fractional derivative can be used to describe the complex problems that involve memory in time because of the nonlocal property.From this, () indicates the memory range at time , and this memory property varies with time.
Using the Caputo fractional derivative, we introduce a memory formalism for two-phase incompressible and immiscible fluid flow in porous media.Denote the wetting phase by a subscript  and the nonwetting phase by .Let   be the saturation of phase .The two-phase saturations are subject to the following constraint: ( For flow in porous media, the velocity u  of each phase  is described by Darcy's law as where K is the absolute permeability tensor in the porous medium,  is the gravity acceleration,  is the depth, and   ,   ,   , and   are the relative permeability, viscosity, pressure, and density of each phase, respectively.The mass conservation equation of each phase is given by where  is the porosity of the medium and   is the external mass flow rate.
The difference between the nonwetting phase and wetting phase pressures is described by the capillary pressure: Denote   =   /  and Φ  =   + .Substituting (3) into (4), we obtain Summing the forms (6) of two phases and taking into account (2), we can reach where   =   +   , Φ  =   + (  −   ) and   =   +   .Furthermore, define u  = −  K∇Φ  and   =   /  .We then have Thus, the equation for the wetting-phase saturation becomes We note that the above two-phase flow formulation is a fractional extension of the formulation [4].We consider the classical mass conservative formulations of incompressible two-phase flow: The porosity often depends on the fluid pressure and saturation.In some problems of incompressible two-phase flow, petroleum reservoir engineering, for example, the porosity varies small with respect to the pressure.However, the injection fluids can change the porosity to a great extent.
In the displacement process of oil by water, the injected water will make the media wet, and then the porosity will be enlarged.Summing the mass conservative formulations of two phases, we obtain Since the porosity varies small with the pressure, we assume / ≃ 0 in (11), and as a result, a pressure equation is obtained with the same form to (9).The wetting-phase saturation equation can be expressed as Because of porosity depending on   , we get by the chain rule that where Υ(  ) =  +     (  ).In general, Υ(  ) can be described by the power law.With this argument, we can deduce the fractional generalization of the classical twophase model.This clearly explains the physical reason to use the fractional model.Finally, we complete our model by the boundary and initial conditions.We divide the boundary Ω of the computational domain Ω into the two nonoverlapping parts: the Dirichelt part Γ  and Newmann part Γ  , where Ω = Γ  ∪ Γ  .The pressure equation ( 6) is subject to the following boundary conditions: where   is the pressure on Γ  , n is the outward unit normal vector to Ω, and u  is the imposed inflow rate on Γ  .The boundary conditions for the saturations are given by The initial saturation of the wetting phase is given by 2.2.Discrete-Fracture Model with Fractional Time Derivatives.Here, the discrete-fracture model [4] is extended to the case with the fractional time derivative.The discretefracture model treats the matrix and fracture gridcells by different geometrical dimensions; that is, if the domain is dimensional, the matrix regions are n-dimensional, but the fractures are simplified as the matrix gridcell interfaces that are ( − 1)-dimensional.This treatment removes the lengthscale contrast resulting from the explicit representation of the fracture aperture as in the single-porosity model, so it is capable to considerably improve the computational efficiency and is convenient in practical implementation.We now decompose the entire domain into two parts: the matrix Ω  and fracture Ω  .The fractures are surrounded by the matrix blocks.We use the subscript  to represent the matrix and the subscript  to represent the fracture system.The pressure in the matrix domain is determined by which is subject to the matrix-fracture interface condition: In the fracture system, we denote the fracture width by  and assume that the potentials are constant along the fracture width, and then obtain the pressure equation in the fracture as where  , is the mass transfer across the matrix-fracture interfaces.The above formulations are similar to the classical model.
As stated previously, the fractional property represents the variability of porosity and permeability, which have a remarkable contrast between matrixes and fractures.As a result, different fractional time derivatives should be used for the flow in the matrix regions and the fracture system.By using the fractional time derivative, we can express the saturation equation in the matrix regions as along with the interface condition Similarly, the saturation equation in the fracture system is given by where  , represents the mass transfer across the matrixfracture interfaces.

Numerical Methods
In this section, we will present the numerical methods for the fractional model of two-phase incompressible flow.In the following, we focus on the time discretization schemes.We firstly divide the total time interval [0, ] into   equal time steps as 0 =  0 <  1 < ⋅ ⋅ ⋅ <    =  and denote the time step length ℎ  = /  .This time division is used for the pressures.Since the saturation varies more rapidly than the pressure, we use a smaller time step size for saturation.Each subinterval (  ,  +1 ] is partitioned into   sub-subintervals as (  ,  +1 ] = ⋃   −1 =0 ( , ,  ,+1 ], where  ,0 =   and  ,  =  +1 , and denote the sub-subinterval length by ℎ  = ( +1 −   )/  .Denote the value of a variable V on the   time point by V  and the one on  , by V , .
For the pressure equation, the saturations take the values of previous time steps, and the capillary potential Φ  on each cell are explicitly calculated by using the cell saturations from the previous time step and the capillary pressure functions.The variables   ,   , and   in the pressure equation are also explicitly calculated by using the cell saturations from Advances in Mathematical Physics the previous time step.From this, we obtain the pressure equation in the matrix domain: where the superscript  represents the time step.Equation ( 17) is subject to the matrix-fracture interface condition: It is similar to express the form in the fracture (referred to by the subscript ) as Once the pressures Φ +1 , and Φ +1 , are computed, the velocities can be evaluated as As previously mentioned, the lower limit  of integral in (1) is a function with time.This function can be chosen according to practical problems.For two-phase flow in porous media, the pressure changes less rapidly than the saturation with the time [18], and hence a large time step is taken for the pressure.We can also see that the pressures are determined by saturations but have not any historical memory.As a result, the memory of saturations can be restricted within each time step for the pressures.For describing this memory property, we define Thus, the Caputo fractional derivative of saturation  is defined by We now introduce the explicit time discretization scheme for approximating where The explicit scheme is employed for the saturation equation both in the matrix domain and in the fracture network For spatial discretization schemes, the cell-centered finite difference method is used for the pressure equation, while the upwind finite volume method is employed for the saturation equation.For the detailed descriptions about spatial discretization schemes, we refer to [23].

Numerical Tests
In this section, two numerical examples are provided to demonstrate the proposed numerical model for two-phase flow with fractional time derivatives.
In all tests, the absolute permeability is a diagonal tensor and the porous media are isotropic.We use the following capillary pressure function [24]: where   is a positive parameter related to the absolute permeability.The relative permeabilities of two phases are computed by We consider a horizontal porous medium of 20 m × 15 m × 1 m with multiple interconnected fractures [23], which is shown in Figure 1.The width of fractures is 0.01 m.The porosities of matrix and fracture media are 0.15 and 1, respectively.The permeabilities in the matrix blocks and the fractures are 50 md and 10 5 md, respectively.The viscosities of the water and oil are all equal to 1 cP.The injection rate is 0.2 PV/year.
Because the medium is horizontal, it is reasonable to neglect the effect of gravity.We inject the water at the left end of the medium, whose void is initially fully saturated with oil, to produce the oil at the right-hand side.There is no other injection and no extraction to the interior of the domain.The fluxes towards outsides of the other boundaries vanish.
Figure 2 shows the effects of fractional time generalization on the average water saturations at different points of PVI.In the legend of Figure 2, the case [  ,   ] represents the orders of the fractional time derivative of water saturation in fracture and matrix regions.From Figure 2, we observe the slower temporal decay when compared to the ordinary case.

Figures 3 , 4 , 5 , 6 , 7 , 8 , 9 ,
10, and 11  show the water saturation contours at different time with three pairs of the fractional time derivative.From these figures, we can see the presence of a time delay effect when the fractional order becomes less than one.This indicates that the porosity and permeability change because of the injected water wetting the media, which makes the fluid flow slowly.
Water saturation contour at 5 years with   = 1 and   = 1.