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

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.

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Γ1−𝛾^𝑡_𝑎(𝑡−𝜏)^−𝛾𝑦′(𝜏)𝑑𝜏, (1) 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: 𝑆_𝑤+𝑆_𝑛=1. (2)

For flow in porous media, the velocity 𝐮_𝛼 of each phase 𝛼 is described by Darcy's law as 𝐮_𝛼𝑘=−^𝑟𝛼𝜇_𝛼𝐊∇𝑝_𝛼+𝜌_𝛼𝑔∇𝑧,𝛼=𝑤,𝑛, (3) where 𝐊 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 𝜙𝐷^𝛾_𝑎,𝑡𝑆_𝛼+∇⋅𝐮_𝛼=𝑞_𝛼,𝛼=𝑤,𝑛, (4) 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: 𝑝_𝑐𝑆_𝑤=𝑝_𝑛−𝑝_𝑤. (5)

Denote 𝜆_𝛼=𝑘_𝑟𝛼/𝜇_𝛼 and Φ_𝛼=𝑝_𝛼+𝜌𝑔𝑧. Substituting (3) into (4), we obtain 𝜙𝐷^𝛾_𝑎,𝑡𝑆_𝛼+∇⋅𝜆_𝛼𝐊∇Φ_𝛼=𝑞_𝛼,𝛼=𝑤,𝑛, (6) Summing the forms (6) of two phases and taking into account (2), we can reach −∇⋅𝜆_𝑡𝐊∇Φ_𝑤−∇⋅𝜆_𝑛𝐊∇Φ_𝑐=𝑞_𝑡, (7) where 𝜆_𝑡=𝜆_𝑤+𝜆_𝑛, Φ_𝑐=𝑝_𝑐+(𝜌_𝑛−𝜌_𝑤)𝑔𝑧 and 𝑞_𝑡=𝑞_𝑤+𝑞_𝑛. Furthermore, define 𝐮_𝑎=−𝜆_𝑡𝐊∇Φ_𝑤 and 𝑓_𝑤=𝜆_𝑤/𝜆_𝑡. We then have 𝐮_𝑤=𝑓_𝑤𝐮_𝑎. (8) Thus, the equation for the wetting-phase saturation becomes 𝜙𝐷^𝛾_𝑎,𝑡𝑆_𝑤+∇⋅𝑓_𝑤𝐮_𝑎=𝑞_𝑤. (9)

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) 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) 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 𝜕𝜙𝑆_𝑤𝑓𝜕𝑡+∇⋅_𝑤𝐮_𝑎=𝑞_𝑤. (12) Because of porosity depending on 𝑆_𝑤, we get by the chain rule that 𝜕𝜙𝑆_𝑤𝑆𝜕𝑡=Υ_𝑤𝜕𝑆_𝑤,𝜕𝑡 (13) where Υ(𝑆_𝑤)=𝜙+𝑆_𝑤𝜙^′(𝑆_𝑤). In general, Υ(𝑆_𝑤) 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 Γ^𝐷 and Newmann part Γ^𝑁, where 𝜕Ω=Γ^𝐷∪Γ^𝑁. The pressure equation (6) is subject to the following boundary conditions: 𝑝_𝑤or𝑝_𝑛=𝑝^𝐷onΓ^𝐷,𝐮_𝛼⋅𝐧=𝐮^𝑁_𝛼onΓ^𝑁, (14) where 𝑝^𝐷 is the pressure on Γ^𝐷, 𝐧 is the outward unit normal vector to 𝜕Ω, and 𝐮^𝑁 is the imposed inflow rate on Γ^𝑁. The boundary conditions for the saturations are given by 𝑆_𝑤=𝑆^𝑁onΓ^𝑁. (15) The initial saturation of the wetting phase is given by 𝑆_𝑤=𝑆⁰_𝑤inΩ. (16)

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 𝑛-dimensional, the matrix regions are 𝑛-dimensional, but the fractures are simplified as the matrix gridcell interfaces that are (𝑛−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 Ω_𝑚 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 −∇⋅𝜆_𝑡,𝑚𝐊_𝑚∇Φ_𝑤,𝑚−∇⋅𝜆_𝑛,𝑚𝐊_𝑚∇Φ_𝑐,𝑚=𝑞_𝑡,𝑚, (17) which is subject to the matrix-fracture interface condition: Φ_𝑤,𝑚=Φ_𝑤,𝑓,Φ_𝑐,𝑚=Φ_𝑐,𝑓,on𝜕Ω_𝑚∩Ω_𝑓. (18) 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) 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 𝜙_𝑚𝐷^𝛾𝑚_𝑎,𝑡𝑆_𝑤,𝑚+∇⋅𝑓_𝑤,𝑚𝐮_𝑎,𝑚=𝑞_𝑤,𝑚, (20) along with the interface condition 𝑆_𝑤,𝑚=𝑆_𝑤,𝑓,on𝜕Ω_𝑚∩Ω_𝑓. (21) Similarly, the saturation equation in the fracture system is given by 𝜙_𝑓𝐷^𝛾𝑓_𝑎,𝑡𝑆_𝑤,𝑓+∇⋅𝑓_𝑤,𝑓𝐮_𝑎,𝑓=𝑞_𝑤,𝑓+𝑄_𝑤,𝑓, (22) where 𝑄_𝑤,𝑓 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,𝑇] into 𝑁_𝑝 equal time steps as 0=𝑡⁰<𝑡¹<⋯<𝑡^𝑁𝑝=𝑇 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 𝑣 on the 𝑡^𝑖 time point by 𝑣^𝑖 and the one on 𝑡^𝑖,𝑗 by 𝑣^𝑖,𝑗.

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 the previous time step. From this, we obtain the pressure equation in the matrix domain: −∇⋅𝜆_𝑡𝑆^𝑖_𝑤,𝑚𝐊_𝑚∇Φ^𝑖+1_𝑤,𝑚−∇⋅𝜆_𝑛𝑆^𝑖_𝑤,𝑚𝐊_𝑚∇Φ^𝑖_𝑐,𝑚=𝑞^𝑖+1_𝑡,𝑚, (23) where the superscript 𝑖 represents the time step. Equation (17) is subject to the matrix-fracture interface condition: Φ^𝑖+1_𝑤,𝑚=Φ^𝑖+1_𝑤,𝑓,Φ^𝑖_𝑐,𝑚=Φ^𝑖_𝑐,𝑓,on𝜕Ω_𝑚∩Ω_𝑓. (24) It is similar to express the form in the fracture (referred to by the subscript 𝑓) as −∇⋅𝜆_𝑡𝑆^𝑖_𝑤,𝑓𝐊_𝑓∇Φ^𝑖+1_𝑤,𝑓−∇⋅𝜆_𝑛𝑆^𝑖_𝑤,𝑓𝐊_𝑓∇Φ^𝑖_𝑐,𝑓=𝑞^𝑖+1_𝑡,𝑓+𝑄^𝑖+1_𝑡,𝑓. (25)

Once the pressures Φ^𝑖+1_𝑤,𝑚 and Φ^𝑖+1_𝑤,𝑓 are computed, the velocities can be evaluated as 𝐮_𝑎^𝑖+1=−𝜆^𝑖_𝑡𝐊∇Φ_𝑤^𝑖+1. (26)

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 𝑎(𝑡)=𝑡^𝑖,𝑡∈𝑡^𝑖,𝑡^𝑖+1. (27) Thus, the Caputo fractional derivative of saturation 𝑆 is defined by 𝐷^𝛾_{𝑡𝑖,𝑡}𝑆(𝑡)=1Γ1−𝛾^𝑡_𝑡^𝑖(𝑡−𝜏)^−𝛾𝜕𝑆(𝜏)𝑡𝜕𝑡𝑑𝜏,𝑡∈^𝑖,𝑡^𝑖+1. (28) We now introduce the explicit time discretization scheme for approximating 𝐷^𝛾_𝑎,𝑡𝑆 at 𝑡^𝑖,𝑗+1 as 𝐷^𝛾_{𝑡𝑖,𝑡𝑖,𝑗+1}𝑆^𝑖,𝑗+1=𝑏₀𝑆^𝑖,𝑗+1−𝑆^𝑖,𝑗𝑏,𝑗=0,₀𝑆^𝑖,𝑗+1−𝑆^𝑖,𝑗+𝑏∑₀^𝑗_𝑘=1𝑏_𝑗−𝑘𝑆^{𝑖,𝑗−𝑘+1}−𝑆^{𝑖,𝑗−𝑘},𝑗=1,…,𝑁_𝑠, (29) where 𝑏₀=ℎ_𝑠^−𝛾Γ,𝑏2−𝛾_𝑘=(𝑘+1)^1−𝛾−𝑘^1−𝛾,𝑘=1,…,𝑁_𝑠−1. (30)

The explicit scheme is employed for the saturation equation both in the matrix domain 𝜙_𝑚𝐷^𝛾𝑚_{𝑡𝑖,𝑡𝑖,𝑗+1}𝑆^𝑖,𝑗+1_𝑤,𝑚+∇⋅𝑓^𝑖,𝑗_𝑤,𝑚𝐮^𝑖+1_𝑎,𝑚=𝑞^𝑖,𝑗+1_𝑤,𝑚, (31) and in the fracture network 𝜙_𝑓𝐷^𝛾𝑓_{𝑡𝑖,𝑡𝑖,𝑗+1}𝑆^𝑖,𝑗+1_𝑤,𝑓+∇⋅𝑓^𝑖,𝑗_𝑤,𝑓𝐮^𝑖+1_𝑎,𝑓=𝑞^𝑖,𝑗+1_𝑤,𝑓+𝑄^𝑖,𝑗_𝑤,𝑓. (32)

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]: 𝑝_𝑐𝑆_𝑤=−𝐵_𝑐log𝑆_𝑤, (33) where 𝐵_𝑐 is a positive parameter related to the absolute permeability. The relative permeabilities of two phases are computed by 𝑘_𝑟𝑤=𝑆³_𝑤,𝑘_𝑟𝑛=1−𝑆_𝑤³. (34)

We consider a horizontal porous medium of 20m×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 10⁵ 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.

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.

Figure 2 

Water saturations with different fractional time derivatives.

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.

Figure 3 

Water saturation contour at 1.25 years with 𝛾𝑓=1 and 𝛾𝑚=1.

Figure 4 

Water saturation contour at 5 years with 𝛾𝑓=1 and 𝛾𝑚=1.

Figure 5 

Water saturation contour at 10 years with 𝛾𝑓=1 and 𝛾𝑚=1.

Figure 6 

Water saturation contour at 1.25 years with 𝛾𝑓=0.9 and 𝛾𝑚=0.95.

Figure 7 

Water saturation contour at 5 years with 𝛾𝑓=0.9 and 𝛾𝑚=0.95.

Figure 8 

Water saturation contour at 10 years with 𝛾𝑓=0.9 and 𝛾𝑚=0.95.

Figure 9 

Water saturation contour at 1.25 years with 𝛾𝑓=0.85 and 𝛾𝑚=0.9.

Figure 10 

Water saturation contour at 5 years with 𝛾𝑓=0.85 and 𝛾𝑚=0.9.

Figure 11 

Water saturation contour at 10 years with 𝛾𝑓=0.85 and 𝛾𝑚=0.9.

Acknowledgments

This work is supported by the Key Program of Shanghai Municipal Education Commission under Grant (no. 12ZZ084) and the Key Project of Chinese Ministry of Education (no. 212109).