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

Authors: Wenwen Zhong; Changpin Li; Jisheng Kou
Journal: Advances in Mathematical Physics
 (2013)

---

## Abstract

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 conventional model does not work well in many cases since it lacks the memory 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.


---

## Body

## 1. 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–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 discrete-fracture 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 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–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.


## 2. 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
(1)Da,tγy(t)=1Γ(1-γ)∫at‍(t-τ)-γy′(τ)dτ,
where y(t) 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 t instead of the constant as usual; that is, a=a(t). The fractional derivative can be used to describe the complex problems that involve memory in time because of the nonlocal property. From this, a(t) indicates the memory range at time t, 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 subscriptw and the nonwetting phase by n. Let Sα be the saturation of phase α. The two-phase saturations are subject to the following constraint:
(2)Sw+Sn=1.For flow in porous media, the velocityuα of each phase α is described by Darcy's law as
(3)uα=-krαμαK(∇pα+ραg∇z),α=w,n,
where K is the absolute permeability tensor in the porous medium, g is the gravity acceleration, z is the depth, and krα,  μα, pα, and ρα are the relative permeability, viscosity, pressure, and density of each phase, respectively.The mass conservation equation of each phase is given by(4)ϕDa,tγSα+∇·uα=qα,α=w,n,
where ϕ is the porosity of the medium and qα is the external mass flow rate.The difference between the nonwetting phase and wetting phase pressures is described by the capillary pressure:(5)pc(Sw)=pn-pw.Denoteλα=krα/μα and Φα=pα+ρgz. Substituting (3) into (4), we obtain
(6)ϕDa,tγSα+∇·λαK∇Φα=qα,α=w,n,
Summing the forms (6) of two phases and taking into account (2), we can reach
(7)-∇·λtK∇Φw-∇·λnK∇Φc=qt,
where λt=λw+λn, Φc=pc+(ρn-ρw)gz and qt=qw+qn. Furthermore, define ua=-λtK∇Φw and fw=λw/λt. We then have
(8)uw=fwua.
Thus, the equation for the wetting-phase saturation becomes
(9)ϕDa,tγSw+∇·(fwua)=qw.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:
(10)∂ϕSα∂t+∇·uα=qα,α=w,n.
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
(11)∂ϕ∂t-∇·λtK∇Φw-∇·λnK∇Φc=qt.
Since the porosity varies small with the pressure, we assume ∂ϕ/∂t≃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
(12)∂ϕSw∂t+∇·(fwua)=qw.
Because of porosity depending on Sw, we get by the chain rule that
(13)∂ϕSw∂t=Υ(Sw)∂Sw∂t,
where Υ(Sw)=ϕ+Swϕ′(Sw). In general, Υ(Sw) can be described by the power law. With this argument, we can deduce the fractional generalization of the classical two-phase 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 ΓD and Newmann part ΓN, where ∂Ω=ΓD∪ΓN. The pressure equation (6) is subject to the following boundary conditions:
(14)pw(orpn)=pDonΓD,uα·n=uαNonΓN,
where pD is the pressure on ΓD, n is the outward unit normal vector to ∂Ω, and uN is the imposed inflow rate on ΓN. The boundary conditions for the saturations are given by
(15)Sw=SNonΓN.
The initial saturation of the wetting phase is given by
(16)Sw=Sw0inΩ.


### 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 discrete-fracture model treats the matrix and fracture gridcells by different geometrical dimensions; that is, if the domain is n-dimensional, the matrix regions are n-dimensional, but the fractures are simplified as the matrix gridcell interfaces that are (n-1)-dimensional. This treatment removes the length-scale 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Ωm and fracture Ωf. The fractures are surrounded by the matrix blocks. We use the subscript m to represent the matrix and the subscript f to represent the fracture system. The pressure in the matrix domain is determined by
(17)-∇·λt,mKm∇Φw,m-∇·λn,mKm∇Φc,m=qt,m,
which is subject to the matrix-fracture interface condition:
(18)Φw,m=Φw,f,Φc,m=Φc,f,on∂Ωm∩Ωf.
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
(19)-∇·λt,fKf∇Φw,f-∇·λn,fKf∇Φc,f=qt,f+Qt,f,
where Qt,f 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(20)ϕmDa,tγmSw,m+∇·(fw,mua,m)=qw,m,
along with the interface condition
(21)Sw,m=Sw,f,on∂Ωm∩Ωf.
Similarly, the saturation equation in the fracture system is given by
(22)ϕfDa,tγfSw,f+∇·(fw,fua,f)=qw,f+Qw,f,
where Qw,f represents the mass transfer across the matrix-fracture interfaces.


## 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
(1)Da,tγy(t)=1Γ(1-γ)∫at‍(t-τ)-γy′(τ)dτ,
where y(t) 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 t instead of the constant as usual; that is, a=a(t). The fractional derivative can be used to describe the complex problems that involve memory in time because of the nonlocal property. From this, a(t) indicates the memory range at time t, 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 subscriptw and the nonwetting phase by n. Let Sα be the saturation of phase α. The two-phase saturations are subject to the following constraint:
(2)Sw+Sn=1.For flow in porous media, the velocityuα of each phase α is described by Darcy's law as
(3)uα=-krαμαK(∇pα+ραg∇z),α=w,n,
where K is the absolute permeability tensor in the porous medium, g is the gravity acceleration, z is the depth, and krα,  μα, pα, and ρα are the relative permeability, viscosity, pressure, and density of each phase, respectively.The mass conservation equation of each phase is given by(4)ϕDa,tγSα+∇·uα=qα,α=w,n,
where ϕ is the porosity of the medium and qα is the external mass flow rate.The difference between the nonwetting phase and wetting phase pressures is described by the capillary pressure:(5)pc(Sw)=pn-pw.Denoteλα=krα/μα and Φα=pα+ρgz. Substituting (3) into (4), we obtain
(6)ϕDa,tγSα+∇·λαK∇Φα=qα,α=w,n,
Summing the forms (6) of two phases and taking into account (2), we can reach
(7)-∇·λtK∇Φw-∇·λnK∇Φc=qt,
where λt=λw+λn, Φc=pc+(ρn-ρw)gz and qt=qw+qn. Furthermore, define ua=-λtK∇Φw and fw=λw/λt. We then have
(8)uw=fwua.
Thus, the equation for the wetting-phase saturation becomes
(9)ϕDa,tγSw+∇·(fwua)=qw.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:
(10)∂ϕSα∂t+∇·uα=qα,α=w,n.
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
(11)∂ϕ∂t-∇·λtK∇Φw-∇·λnK∇Φc=qt.
Since the porosity varies small with the pressure, we assume ∂ϕ/∂t≃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
(12)∂ϕSw∂t+∇·(fwua)=qw.
Because of porosity depending on Sw, we get by the chain rule that
(13)∂ϕSw∂t=Υ(Sw)∂Sw∂t,
where Υ(Sw)=ϕ+Swϕ′(Sw). In general, Υ(Sw) can be described by the power law. With this argument, we can deduce the fractional generalization of the classical two-phase 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 ΓD and Newmann part ΓN, where ∂Ω=ΓD∪ΓN. The pressure equation (6) is subject to the following boundary conditions:
(14)pw(orpn)=pDonΓD,uα·n=uαNonΓN,
where pD is the pressure on ΓD, n is the outward unit normal vector to ∂Ω, and uN is the imposed inflow rate on ΓN. The boundary conditions for the saturations are given by
(15)Sw=SNonΓN.
The initial saturation of the wetting phase is given by
(16)Sw=Sw0inΩ.


## 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 discrete-fracture model treats the matrix and fracture gridcells by different geometrical dimensions; that is, if the domain is n-dimensional, the matrix regions are n-dimensional, but the fractures are simplified as the matrix gridcell interfaces that are (n-1)-dimensional. This treatment removes the length-scale 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Ωm and fracture Ωf. The fractures are surrounded by the matrix blocks. We use the subscript m to represent the matrix and the subscript f to represent the fracture system. The pressure in the matrix domain is determined by
(17)-∇·λt,mKm∇Φw,m-∇·λn,mKm∇Φc,m=qt,m,
which is subject to the matrix-fracture interface condition:
(18)Φw,m=Φw,f,Φc,m=Φc,f,on∂Ωm∩Ωf.
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
(19)-∇·λt,fKf∇Φw,f-∇·λn,fKf∇Φc,f=qt,f+Qt,f,
where Qt,f 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(20)ϕmDa,tγmSw,m+∇·(fw,mua,m)=qw,m,
along with the interface condition
(21)Sw,m=Sw,f,on∂Ωm∩Ωf.
Similarly, the saturation equation in the fracture system is given by
(22)ϕfDa,tγfSw,f+∇·(fw,fua,f)=qw,f+Qw,f,
where Qw,f represents the mass transfer across the matrix-fracture interfaces.


## 3. 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,T] into Np equal time steps as 0=t0<t1<⋯<tNp=T and denote the time step length hp=T/Np. 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 (ti,ti+1] is partitioned into Ns sub-subintervals as (ti,ti+1]=⋃j=0Ns-1‍(ti,j,ti,j+1], where ti,0=ti and ti,Ns=ti+1, and denote the sub-subinterval length by hs=(ti+1-ti)/Ns. Denote the value of a variable v on the ti time point by vi and the one on ti,j by vi,j.For the pressure equation, the saturations take the values of previous time steps, and the capillary potentialΦc on each cell are explicitly calculated by using the cell saturations from the previous time step and the capillary pressure functions. The variables λw, λn, and λt in the pressure equation are also explicitly calculated by using the cell saturations from the previous time step. From this, we obtain the pressure equation in the matrix domain:
(23)-∇·λt(Sw,mi)Km∇Φw,mi+1-∇·λn(Sw,mi)Km∇Φc,mi=qt,mi+1,
where the superscript i represents the time step. Equation (17) is subject to the matrix-fracture interface condition:
(24)Φw,mi+1=Φw,fi+1,Φc,mi=Φc,fi,on∂Ωm∩Ωf.
It is similar to express the form in the fracture (referred to by the subscript f) as
(25)-∇·λt(Sw,fi)Kf∇Φw,fi+1-∇·λn(Sw,fi)Kf∇Φc,fi=qt,fi+1+Qt,fi+1.Once the pressuresΦw,mi+1 and Φw,fi+1 are computed, the velocities can be evaluated as
(26)uai+1=-λtiK∇Φwi+1.As previously mentioned, the lower limita 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
(27)a(t)=ti,t∈(ti,ti+1].
Thus, the Caputo fractional derivative of saturation S is defined by
(28)Dti,tγS(t)=1Γ(1-γ)∫tit‍(t-τ)-γ∂S∂t(τ)dτ,t∈(ti,ti+1].
We now introduce the explicit time discretization scheme for approximating Da,tγS at ti,j+1 as
(29)D~ti,ti,j+1γSi,j+1={b0(Si,j+1-Si,j),j=0,b0(Si,j+1-Si,j)+b0∑k=1j‍bj-k(Si,j-k+1-Si,j-k),j=1,…,Ns,
where
(30)b0=hs-γΓ(2-γ),bk=(k+1)1-γ-k1-γ,k=1,…,Ns-1.The explicit scheme is employed for the saturation equation both in the matrix domain(31)ϕmD~ti,ti,j+1γmSw,mi,j+1+∇·(fw,mi,jua,mi+1)=qw,mi,j+1,
and in the fracture network
(32)ϕfD~ti,ti,j+1γfSw,fi,j+1+∇·(fw,fi,jua,fi+1)=qw,fi,j+1+Qw,fi,j.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].


## 4. 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]:
(33)pc(Sw)=-Bclog(Sw),
where Bc is a positive parameter related to the absolute permeability. The relative permeabilities of two phases are computed by
(34)krw=Sw3,krn=(1-Sw)3.We consider a horizontal porous medium of20m×15m×1m 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 105 md, respectively. The viscosities of the water and oil are all equal to 1 cP. The injection rate is 0.2 PV/year.Figure 1
Distribution of fractures.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.Figure2 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 [γf,γm] 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.Figure 2
Water saturations with different fractional time derivatives.Figures3, 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.Figure 3
Water saturation contour at 1.25 years withγf=1 and γm=1.Figure 4
Water saturation contour at 5 years withγf=1 and γm=1.Figure 5
Water saturation contour at 10 years withγf=1 and γm=1.Figure 6
Water saturation contour at 1.25 years withγf=0.9 and γm=0.95.Figure 7
Water saturation contour at 5 years withγf=0.9 and γm=0.95.Figure 8
Water saturation contour at 10 years withγf=0.9 and γm=0.95.Figure 9
Water saturation contour at 1.25 years withγf=0.85 and γm=0.9.Figure 10
Water saturation contour at 5 years withγf=0.85 and γm=0.9.Figure 11
Water saturation contour at 10 years withγf=0.85 and γm=0.9.



---

*Source: 429835-2013-08-06.xml*