Advances in Mathematical Physics

Volume 2013 (2013), Article ID 429835, 7 pages

http://dx.doi.org/10.1155/2013/429835

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

^{1}Department of Mathematics, Shanghai University, Shanghai 200444, China^{2}School of Mathematics and Statistics, Hubei Engineering University, Xiaogan, Hubei 432000, China

Received 17 May 2013; Revised 13 July 2013; Accepted 13 July 2013

Academic Editor: H. Srivastava

Copyright © 2013 Wenwen Zhong et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 where is a function and . 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 of each phase is described by Darcy's law as 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 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 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 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 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: 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 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 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 ()-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 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 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 into equal time steps as 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 is partitioned into sub-subintervals as , where and , and denote the sub-subinterval length by . 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: 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 and 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 at as 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].

#### 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]: 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 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 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.

#### 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).

#### References

- R. G. Baca, R. C. Arnett, and D. W. Langford, “Modelling fluid flow in fractured-porous rock masses by finite-element techniques,”
*International Journal for Numerical Methods in Fluids*, vol. 4, no. 4, pp. 337–348, 1984. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - G. I. Barenblatt, I. P. Zheltov, and I. N. Kochina, “Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata],”
*Journal of Applied Mathematics and Mechanics*, vol. 24, no. 5, pp. 1286–1303, 1960. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - Z. Chen, G. Huan, and Y. Ma,
*Computational Methods for Multiphase Flows in Porous Media*, SIAM, Philadelphia, Pa, USA, 2006. View at Publisher · View at Google Scholar · View at MathSciNet - H. Hoteit and A. Firoozabadi, “An efficient numerical model for incompressible two-phase flow in fractured media,”
*Advances in Water Resources*, vol. 31, no. 6, pp. 891–905, 2008. View at Publisher · View at Google Scholar · View at Scopus - H. Kazemi, “Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution,”
*SPE Journal*, vol. 9, no. 4, pp. 451–462, 1969. View at Publisher · View at Google Scholar · View at Scopus - S. H. Lee, C. L. Jensen, and M. F. Lough, “Efficient finite-difference model for flow in a reservoir with multiple length-scale fractures,”
*SPE Journal*, vol. 5, no. 3, pp. 268–275, 2000. View at Google Scholar · View at Scopus - J. Noorishad and M. Mehran, “An upstream finite element method for solution of transient transport equation in fractured porous media,”
*Water Resources Research*, vol. 18, no. 3, pp. 588–596, 1982. View at Publisher · View at Google Scholar · View at Scopus - K. Pruess and T. N. Narasimhan, “A practical method for modeling fluid and heat flow in fractured porous media,”
*SPE Journal*, vol. 25, no. 1, pp. 14–26, 1985. View at Google Scholar · View at Scopus - S. Sun, A. Firoozabadi, and J. Kou, “Numerical modeling of two-phase binary fluid mixing using mixed finite elements,”
*Computational Geosciences*, vol. 16, no. 4, pp. 1101–1124, 2012. View at Publisher · View at Google Scholar - L. K. Thomas, T. N. Dixon, and R. G. Pierson, “Fractured reservoir simulation,”
*SPE Journal*, vol. 23, no. 1, pp. 42–54, 1983. View at Publisher · View at Google Scholar · View at Scopus - J. E. Warren and P. J. Root, “The behavior of naturally fractured reservoirs,”
*SPE Journal*, vol. 3, no. 3, pp. 245–255, 1963. View at Publisher · View at Google Scholar - E. di Giuseppe, M. Moroni, and M. Caputo, “Flux in porous media with memory: models and experiments,”
*Transport in Porous Media*, vol. 83, no. 3, pp. 479–500, 2010. View at Publisher · View at Google Scholar · View at Scopus - M. Caputo, “Models of flux in porous media with memory,”
*Water Resources Research*, vol. 36, no. 3, pp. 693–705, 2000. View at Publisher · View at Google Scholar · View at Scopus - R. Garra, “Fractional-calculus model for temperature and pressure waves in fluid-saturated porous rocks,”
*Physical Review E*, vol. 84, no. 3, Article ID 036605, 6 pages, 2011. View at Publisher · View at Google Scholar - K. Diethelm, N. J. Ford, A. D. Freed, and Yu. Luchko, “Algorithms for the fractional calculus: a selection of numerical methods,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 194, no. 6–8, pp. 743–773, 2005. View at Publisher · View at Google Scholar · View at Scopus - H. Ding and C. Li, “Mixed spline function method for reaction-subdiffusion equations,”
*Journal of Computational Physics*, vol. 242, pp. 103–123, 2013. View at Publisher · View at Google Scholar - C. Li, A. Chen, and J. Ye, “Numerical approaches to fractional calculus and fractional ordinary differential equation,”
*Journal of Computational Physics*, vol. 230, no. 9, pp. 3352–3368, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - Z. Chen, G. Huan, and B. Li, “An improved IMPES method for two-phase flow in porous media,”
*Transport in Porous Media*, vol. 54, no. 3, pp. 361–376, 2004. View at Publisher · View at Google Scholar · View at MathSciNet - X.-J. Yang, H. M. Srivastava, J.-H. He, and D. Baleanu, “Cantor-type cylindrical-coordinate method for differential equations with local fractional derivatives,”
*Physics Letters A*, vol. 377, no. 28–30, pp. 1696–1700, 2013. View at Publisher · View at Google Scholar - K. H. Coats, “Reservoir simulation: state-of-the-art,”
*Journal of Petroleum Technology*, vol. 34, no. 8, pp. 1633–1642, 1982. View at Publisher · View at Google Scholar - R. G. Fagin and C. H. Stewart Jr., “A new approach to the two-dimensional multiphase reservoir simulator,”
*SPE Journal*, vol. 6, no. 2, pp. 175–182, 1966. View at Publisher · View at Google Scholar - P. A. Forsyth Jr. and P. H. Sammon, “Quadratic convergence for cell-centered grids,”
*Applied Numerical Mathematics*, vol. 4, no. 5, pp. 377–394, 1988. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - J. Kou, S. Sun, and B. Yu, “Multiscale time-splitting strategy for multiscale multiphysics processes of two-phase flow in fractured media,”
*Journal of Applied Mathematics*, vol. 2011, Article ID 861905, 24 pages, 2011. View at Publisher · View at Google Scholar · View at MathSciNet - H. Hoteit and A. Firoozabadi, “Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures,”
*Advances in Water Resources*, vol. 31, no. 1, pp. 56–73, 2008. View at Publisher · View at Google Scholar · View at Scopus