Research Article  Open Access
Jisheng Kou, Shuyu Sun, Bo Yu, "Multiscale TimeSplitting Strategy for Multiscale Multiphysics Processes of TwoPhase Flow in Fractured Media", Journal of Applied Mathematics, vol. 2011, Article ID 861905, 24 pages, 2011. https://doi.org/10.1155/2011/861905
Multiscale TimeSplitting Strategy for Multiscale Multiphysics Processes of TwoPhase Flow in Fractured Media
Abstract
The temporal discretization scheme is one important ingredient of efficient simulator for twophase flow in the fractured porous media. The application of singlescale temporal scheme is restricted by the rapid changes of the pressure and saturation in the fractured system with capillarity. In this paper, we propose a multiscale time splitting strategy to simulate multiscale multiphysics processes of twophase flow in fractured porous media. We use the multiscale time schemes for both the pressure and saturation equations; that is, a large timestep size is employed for the matrix domain, along with a small timestep size being applied in the fractures. The total time interval is partitioned into four temporal levels: the first level is used for the pressure in the entire domain, the second level matching rapid changes of the pressure in the fractures, the third level treating the response gap between the pressure and the saturation, and the fourth level applied for the saturation in the fractures. This method can reduce the computational cost arisen from the implicit solution of the pressure equation. Numerical examples are provided to demonstrate the efficiency of the proposed method.
1. Introduction
The fractured porous media is composed of two distinct pore spaces: the fracture and the matrix. The fracture spaces have higher permeability than the matrix, but possess very small volume compared with the matrix. Numerical simulation of multiphase flow in fractured porous media is of importance in the subsurface flow, environmental problems and petroleum reservoir engineering. In this paper, we focus on twophase fluid flow in fractured porous media with capillarity. In the displacement of the nonwetting phase by the wetting phase, the injected fluid flows rapidly across the fracture network if the effect of capillarity is neglected; however, when the capillarity is taken into account, the capillary contrast between the matrix and the fracture may have a significant influence on flow paths. In the fractured porous media, the different capillary pressure functions are employed within the fracture and the matrix blocks because of their different permeability types; that is, the capillary pressure functions are not continuous on the matrixfracture interfaces. Therefore, the discontinuity of the saturation is predicted because of the continuity of capillary pressure [1].
To simulate the multiphase flow in fractured porous media, several conceptual models have been developed [2–16], such as the singleporosity model, the dualporosity/dualpermeability model, and the discretefracture model. In the singleporosity model, the fractures are treated similarly to the matrix. This model requires the excessive gridcells because of the contrast of geometrical scales between the matrix and the fracture. The dualporosity/dualpermeability model [3, 8–10, 13, 15, 16] envisions the matrix blocks and the fractures as two separate parts, which are interconnected by fluid mass transfer across their interfaces. The discrete fracture model [2, 11, 12] is based on the conception that the fracture aperture is very small compared to the matrix blocks, and as a result, we can simplify the fracture as the lower dimensional domain to reduce the contrast of geometric scales occurring in the singleporosity model. The discrete fracture model is shown to be is numerically superior to the singleporosity model and the dualporosity models [6, 7]. This model will be used in this paper.
The implicit pressure explicit saturation (IMPES) method [17–24] is a popular timestepping approach employed in multiphase flow simulation. Based on the property of multiphysics processes, the IMPES method splits the coupled system into one pressure equation and one saturation equation. The saturation and capillary pressure in the pressure equation are treated explicitly to eliminate its nonlinearity. These treatments result in the decoupling between the pressure equation and the saturation equation, and the instability of the IMPES method [21], especially for the case with capillarity in the heterogeneous media. In order to improve the stability of IMPES, we present an implicit treatment of capillarity in [25], which will be described in Section 3. The iterative version of this method is presented in [26]. Based on the previous works, we will discuss this approach being applied in the fractured media. For the fluid flow in porous media, the pressure changes less rapidly than the saturation with the time, and hence we can take a larger timestep for the pressure than the one for saturation. This approach is an improvement of IMPES originally proposed in [17, 24, 27]. It is also demonstrated that this technique is very efficient for the method proposed in [25]. Relevant work in coupled flow and transport simulation can be found in, for example, [28–35]. In this paper, we will describe a multiscale time strategy for multiphysics processes of twophase flow in fractured media.
For fractured media, the pressure and saturation in the fracture network change more rapidly than the ones in the matrix blocks. The explicit treatment of saturation in the pressure equation may be instable, and thus the small timestep size is required if the same time scale is employed in the entire spatial domain. The explicit scheme used for updating the saturation suffers from Courant number stability constraints, and as a result, we should use the smallest timestep size to match the rapid changes in the fracture network if using the singlescale timestepping. To resolve this problem, we should use the multiscale timestepping for the pressure and the saturation in the fracture and the matrix.
Multiscale time splitting strategy has been studied in the literature, for example, [36–45]. An explicit subtiming scheme is proposed in [41, 42]. The computational domain is divided into separate coarse and finegrid locations. We fist carry out the computations for all coarse gridblocks by using a timestep size that satisfies the coarse gridblock stability constraint. Once the solutions in a coarsely discretized region are obtained, the unknowns within the fine grids are solved by using a smaller timestep size that matches the fine grid stability constraints. In the subtiming methodology for implicittype timestepping schemes proposed in [39], the fully implicit schemes with multiple timestep sizes are utilized to the different portions of a domain, which have different accuracy requirements. This technique has been applied to the timedependent flow and transport simulations in fractured porous media [40], which uses the implicit subtimestepping in the locations within and near the fractures. It is demonstrated that the implicit multiple timestepping is capable to significantly improve the accuracy and efficiency of numerical simulation. The hybrid implicitexplicit approach [46] treats implicitly the portions of the domain with the particular stability requirements, along with the explicit timestepping applied to the remainder of the domain. The advantage of this approach is that it can reduce the size of the problems being solved implicitly at each timestep.
In this paper, we present a multiscale time splitting strategy for twophase flow in fractured porous media. Taking account into the different physical properties of the matrix blocks and the fracture network, and the multiphysics processes of twophase flow, we apply multiple time scales for the pressure and saturation equations. The pressure equation is formed from the conception proposed in [25]; that is, the capillary pressure is treated implicitly in the pressure equation. For the pressure, we employ a large timestep size in the matrix domain, along with a small timestep size being applied in the fracture system. The timestep for the pressure in the fractures is further partitioned to a few subtimesteps, in which the saturation will be updated explicitly. The updating process of the saturation employs twoscale timestepping; a coarse timestepping is applied in the matrix blocks, along with a fine one in the fracture system.
The rest of this paper is organized as follows. In Section 2, we review the twophase incompressible flow model, along with the discretefracture model. In Section 3, we describe our multiscale time splitting strategy. Here, the cellcentered finite difference (CCFD) method [47] is employed for spatial discretization, but it can be extended to the other discretization schemes. In Section 4, some numerical examples are given to demonstrate the advantage of the proposed method. Finally, we summarize this work.
2. Mathematical Model of TwoPhase Incompressible Flow
2.1. Governing Equations of TwoPhase Flow
We now describe twophase incompressible and immiscible fluid flow in porous media, and denote the wetting phase by a subscript and the nonwetting phase by . The mass conservation equation of each phase is given by where is the porosity of the medium, is the saturation of phase , and is the external mass flow rate. In (2.1), Darcy's velocity of each phase is given by 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.
In addition, the saturations of two phases are subject to the following constraint The difference between the nonwetting phase and wetting phase pressures is described by the capillary pressure
In this paper, we use the twophase flow formulation [1, 7, 48] given by where , , , , , , , and . From the above definitions, the wettingphase velocity is expressed as
Finally, we consider 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 (2.5) 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. DiscreteFracture Model
Here, we employ the discretefracture model [6, 7] to explicitly describe the fractures in fractured porous media. In this model, different geometrical dimensions are taken for the matrix and fracture grid cells; that is, for dimensional domain, the matrix gridcells are dimensional, but the fracture gridcells are simplified as the matrix gridcell interfaces that are ()dimensional. This treatment is capable to considerably improve the computational efficiency since it can remove the lengthscale contrast resulting from the explicit representation of the fracture aperture as in the singleporosity model. The entire domain is decomposed into two parts: the matrix (dimensional) and fracture (()dimensional), and the fractures are surrounded by the matrix blocks. The pressure equation in the matrix domain is given by where the subscript represents the matrix domain. Equation (2.11) is subject to the matrixfracture interface condition The fracture width is denoted by . In the fracture system, we assume the potentials and saturations are constant along the fracture width, and then by integrating the equation, obtain the pressure equation in the fracture (represented by the subscript ) as where is the mass transfer across the matrixfracture interfaces. Similarly, we can express the saturation equation in the matrix domain as with the interface condition In accordance, the saturation equation in the fracture system is given by where represents the mass transfer across the matrixfracture interfaces. This form of saturation equation will be used to rebuild the coupling relationship between the pressure and saturation. The other form of saturation equation is used to compute the saturation, which is expressed as in the matrix domain and in the fractures where is the same to , but with different expressions in the spatial discretization.
3. MultiScale Time Splitting Strategy
In this paper, we use the cellcentered finite difference (CCFD) method for the spatial discretization, but apparently, our approach can be also extended to the other spatial discretization schemes. In this paper, we focus on the case of twodimensional space. The extension to threedimensional case is straightforward. We assume that the porous media are isotropic and the absolute permeability is a scalar function.
3.1. MultiScale Time Splitting Approach for the Pressure Equation with Implicit Capillary Pressure
We firstly divide the total time interval into timesteps as . This time division is used for the pressure in the matrix domain. Define the timestep length . Since the pressure in the fractures varies more rapidly than the matrix domain, we use a smaller timestep size for the pressure in the fractures; that is, we partition each subinterval into subsubintervals as , where and , and denote the subsubinterval length by . We denote the value of a variable on the time point by and the one on by .
We consider the grid system with rectangular cells for the matrix domain. The fracture that is onedimensional can be partitioned into many spatial subintervals. The classical IMPES method is instable because each cell capillary potential is explicitly calculated by using the cell saturations from the previous timestep and the capillary pressure functions. In this paper, we employ the implicit approach to handle capillary pressure proposed in [25]. This method treats the variables , , and in the pressure equation as IMPES, but the capillary potentials will be treated implicitly. From this, we obtain the mass conservation equation in the matrix domain where the superscript represents the timestep and is the temporal discretization parameter. Apparently, it becomes the form of IMPES if . Here, we always take ; that is, the capillarity is implicitly treated. It is similar to obtain the form in the fracture (referred to by the subscript ) as
Let indicate the boundary of the matrix gridcells. We use the CCFD scheme applied to (3.1) and obtain where indicates the area of the cell , and indicates the flux across the boundary of the gridcell . Here, we assume that is also a gridcell of the fractures if . On each interface , is the unit normal vector pointing from to . If and , the fluxes in (3.3) are given by where is the length of and stands for the Euclidean distance between the central points of the cells and . In (3.4), the terms and are given by where stands for the Euclidean distance from the central points of the cell and the cell boundary . If and is a gridcell of the fracture system, we have where and are given by For the case locating on the entire domain boundary, the pressure potential is provided by Dirichelt boundary conditions if , otherwise the fluxes are determined by the Newmann conditions if .
The above processes of deriving (3.3)–(3.9) can also be carried out for the reduceddimensional fracture system. Let be a gridcell of the fracture network. The CCFD scheme is applied to (3.2), and then we have where indicates the face of the gridcell in the fracture system and indicates the flux across the boundary of the fracture gridcell . Different from (3.3), no potential of the matrix domain is required for computing the fracturefracture fluxes in (3.10), but the matrixfracture transfer needs to be treated as a source term in the fracture system. We consider a gridcell of the fracture network where and and are the gridcells of the matrix domain. The fracture width is denoted by . The volumetric transfer across the matrixfracture interface is given by where , , and are defined by (3.6)–(3.8), respectively. We now consider the treatment of the interface connecting multiple fractures. Let be the set composed of the fracture grid cells that are jointed by . The discretization of the mass conservation equation is given by where and .
Combining the formulations in the matrix domain and fracture network, we obtain the discretization of total mass conservation equation given by The terms , , , and of (3.13) arise from the interconnections of the matrix blocks and fracture system. Since and is unknown, the above equation can not be solved immediately because of its nonlinearity. In order to eliminate the nonlinearity, we approximate the capillary pressure at by Furthermore, we employ the backward Euler time discretization for the saturation equation both in the matrix domain and in the fracture network It is analogous to the derivation of the pressure equation (3.13), and we approximate (3.15) and (3.16) by the CCFD method as Substituting (3.14) and (3.17) into (3.13), we obtain the coupled pressure equation where
Equation (3.18), along with (3.19), is utilized to compute the matrix pressure in a large timestep. In order to match the rapid changes of the pressure in the fracture network, we use a small timestep size in the fracture system. For each time substep, we consider the pressure equation in the fractures as It is obtained by using the CCFD scheme to (3.20) that where is the face of the gridcell in the fracture system and represents the flux across the boundary of the fracture gridcell . Let and are the matrix gridcells and , then in (3.21), the matrixfracture transfer is given by where along with The treatment of the interface connecting multiple fractures is similar to (3.12). Therefore, in the timestep used for the pressure in the fractures, we solve the following pressure equation where
In the above equations, the inverse of and is not expensive as they are diagonal. In the timestep , we solve the linear system (3.18) implicitly to obtain the pressure , and then in the time substep compute by solving (3.26).
Once the pressures and is obtained, the fluxes can be evaluated as described below. For a boundary of matrix gridcell , is the unit normal vector pointing towards outside . If and , where If and is a fracture gridcell, the fluxes are computed by (3.23). The fluxes on the domain boundary are obtain by the boundary conditions as discussed above.
3.2. MultiScale Explicit TimeStepping for the Saturation Equation
As the pressure computed by the abovedescribed method has coupled with the saturation, the forward Euler time discretization is used for computing the saturation equation of the wetting phase in each timestep used for the pressure in the fracture. According to the concept proposed in [17, 24, 27], the timestep size for the pressure can be taken larger than the one for saturation, which is also demonstrated to be efficient for the method proposed in [25]. We divide the timestep for the pressure in the fracture network into timesteps as , where and . This time division is used for the saturation in the matrix domain. Because of the permeability contrast between the matrix and the fracture, the saturation in the fractures changes more rapidly than that in the matrix domain. Therefore, a smaller timestep size is required for the saturation in the fracture system; that is, we partition into time substeps as , where and . Denote and . The explicit scheme is employed for the saturation equation both in the matrix domain and in the fracture network We use the upwind CCFD method to discretize the saturation equation (3.30) Let be the interface between the matrix gridcells and ; that is, . If , the term in (3.32) is determined by If is a gridcell of the fracture network, the term in (3.32) is determined by where
It is analogous to discretize the saturation equation in the fracture system as where is a gridcell of the fracture network. Let where and are the fracture gridcells, then we have The volumetric transfer across the matrixfracture interface is given by where In the above equations, the fluxes are obtained after solving the pressure equation.
To explicitly update of the saturation in the fractures, we use the previous matrix saturation in the matrixfracture interfaces, and then obtain the matrixvector form of (3.36) After smaller timesteps, we update the saturation in the matrix domain by employing the following matrixvector form of (3.32) as
In (3.40) and (3.41), and indicate the interconnections of the saturation in the matrix blocks and fracture system.
4. Numerical Tests
Four numerical examples are provided to demonstrate the efficiency of the multiscale time splitting strategy proposed in this work. We compare it with a conventional method that takes the larger timestep for the pressure than the saturation. The conventional method is referred to the method proposed in [25] when the capillarity is taken into account, while it becomes the classical IMPES method for the case without the effect of capillarity. The two methods use the same relaxation factor for each designated example with the capillarity.
4.1. Physical and Computational Data Used in Numerical Tests
The normalized saturation is given by where and are the residual saturations of the wetting and nonwetting phases, respectively. The capillary pressure function [1] is given by where is a positive parameter related to the absolute permeability. The relative permeabilities of two phases are given by where is a positive integer number.
In Examples 1 and 2, we consider the horizontal layers with the same domain dimensions 10 m × 10 m × 1 m, but containing the different networks composed of two fractures. The domain of Examples 3 and 4 is a horizontal porous medium of 20 m × 15 m × 1 m with multiple interconnected fractures. With media being horizontal, it is reasonable to neglect the effect of gravity in all examples. We simulate the displacement process of oil by water; that is, we inject the water at the left end of the medium, which void is initially fully saturated with oil, to produce the oil at the righthand side. The fluxes towards outsides of the other boundaries vanish. There is no other injection and no extraction to the interior of the domain.
The physical an computational data is provided in Table 1, in which the subscripts "" and "" of some raw titles represent the matrix domain and the fracture system, respectively. In Tables 2–5, the title for column represents the timestep size used for the pressure in the matrix domain if the proposed method is applied, and the other column titles are defined as in Section 3.





4.2. Example 1
In Example 1, the fracture distribution in the tested medium is provided in Figure 1. The permeabilities in the matrix blocks and the fractures are 1 md and 10^{5} md, respectively. The viscosities of the water and oil are 1 cP and 0.5 cP, respectively. The injection rate is 0.2 PV/year. We use the cubic relative permeabilities. The total number of matrix and fracture gridcells is 2704 in this example. The residual saturations of water and oil are zero. The other data used in this example is provided in Table 1.
We simulate the displacement process of oil by water the saturation until 0.35 pore volume injection (PVI). Displayed in Table 2 are the timestep sizes and multiple subtiming strategy required for obtaining the stable solutions by the proposed multiscale time strategy and the conventional temporal approach. In Figures 2 and 3, we show the water saturation profiles at 0.35 PVI computed by the two methods.
In this example, the twoscale time splitting strategy is used only for the pressure equation, but with the singlescale timestepping for the saturation equation. The two methods have the same subtiming between the pressure and the saturation, that is, . In order to achieve stable solutions, with twoscale strategy for the pressure equation, the proposed method can take the timestep size 0.9827 days, which is twice as many over the maximum timestep size 0.4121 days required by the conventional method. It is evident that the size of the pressure equation in the finescale timesteps is less than that in the coarse scale. Hence, this example demonstrates the efficiency of the proposed twoscale implicit treatment for the pressure equation.
4.3. Example 2
In this example, we consider a fractured porous medium with two interconnected fractures, which is shown in Figure 4. The permeability in the fractures is 10^{6} md. The viscosity of the oil is 0.65 cP. In this example, the total number of matrix and fracture gridcells is 2704. The other data used in this example is provided in Table 1. The displacement of oil by water is simulated until 0.45 PVI.
In this example, the proposed method takes the twoscale time splitting strategy for both the pressure and the saturation. The same subtiming between the pressure and the saturation are also employed in the two methods. The details are displayed in Table 3. The water saturation profiles at 0.45 PVI computed by the two methods are shown in Figures 5 and 6, respectively.
As shown in Table 3, in order to achieve stable solutions, the proposed method, with twoscale subtiming, can take a very large timestep size 1.6846 days, while the conventional method requires a very small timestep size, which is less than a tenth part of the proposed method. This example demonstrates that the twoscale time splitting strategy is efficient not only for the pressure equation, but also for the saturation equation.
4.4. Example 3
In this example, we consider a porous medium with a network composed of multiple interconnected fractures, which is shown in Figure 7. In Example 3, the total number of matrix and fracture gridcells is 3300. The commotional data used in this example is provided in Table 1. We simulate the displacement of oil by water until 0.5 PVI.
In this example, we use the proposed method with twoscale time splitting strategy for both the pressure and the saturation, along with the subtiming approach between the pressure and the saturation that is also employed in the conventional method. The details are provided in Table 4. Figures 8 and 9 show the water saturation profiles at 0.5 PVI computed by the two methods, respectively.
From Table 4, we can see that the multiscale time splitting technique is capable of taking a very large timestep size 1.9010 days to achieve stable solutions. The conventional method, however, requires a very small timestep size, that is, 0.0895 days, which is the maximum timestep size for its stability with . Apparently, the multiscale time splitting strategy is superior to the conventional method for the multiplefractured media.
4.5. Example 4
In this example, we attempt to provide the performance of the proposed method for the case neglecting the capillarity. The domain, fluid and computational parameters of this example are the same to Example 3, but the capillarity is absent in this example. The displacement of oil by water is also simulated until 0.5 PVI.
Displayed in Table 5 is the comparison of the proposed multiscale time splitting strategy and the conventional method. Figures 10 and 11 are the water saturation profiles at 0.5 PVI computed by the two methods, respectively. From Figures 10 and 11, we can see that without the effect of capillary pressure contrast at the matrixfracture interfaces, the water flows very quickly through the fractures. The results in Table 5 indicate that the proposed multiscale time splitting method is still efficient.
5. Conclusions
A multiscale time splitting strategy is introduced for simulating twophase flow in fractured porous media. In the fractured porous media, the fracture system is high permeable, but with the very small storage. With different capillary pressure functions being imposed in the matrix blocks and the fracture network, the discontinuity of the saturation on the matrixfracture interface results from the continuity of capillary pressure. As a result, the conventional singlescale temporal schemes have shortage to handle the complexity of the multiphysics processes of twophase flow. In this paper, we use multiscale time schemes for both the pressure and saturation equations. Based on the conception proposed in [25], the capillary pressure is treated implicitly in the pressure equation, but the saturation is yet explicitly updated. Our multiscale time splitting approach divides the total time interval into four temporal levels; the first level is used for the pressure in the entire domain, the second level matches the rapid changes of the pressure in the fractures, the third level is used to treat the response gap between the pressure and saturation, and the fourth level is used to match the fast changes of the saturation in the fracture system.
Numerical (Examples 1–3) demonstrate that when the effect of capillarity is considered, the multiscale time splitting approach can take a very large timestep size to achieve stable solutions, but a very small timestep size is required by the conventional method. The computational efficiency can be improved as it can reduce the cost of implicit solution of the pressure equation. Example 4 indicates that the proposed method is also efficient for the case without capillarity.
Acknowledgment
The authors cheerfully appreciate the generous support of the university research fund to the Computational Transport Phenomena Laboratory at KAUST.
References
 H. Hoteit and A. Firoozabadi, “Numerical modeling of twophase flow in heterogeneous permeable media with different capillarity pressures,” Advances in Water Resources, vol. 31, no. 1, pp. 56–73, 2008. View at: Publisher Site  Google Scholar
 R. G. Baca, R. C. Arnett, and D. W. Langford, “Modelling fluid flow in fracturedporous rock masses by finiteelement techniques,” International Journal for Numerical Methods in Fluids, vol. 4, no. 4, pp. 337–348, 1984. View at: Google Scholar
 G. Barenblatt, Y. Zheltov, and I. Kochina, “Basic concepts in the theory of seepage of homogeneous fluids in fissurized rocks,” Journal of Applied Mathematics and Mechanics, vol. 24, no. 5, pp. 1286–1303, 1983. View at: Google Scholar
 Z. Chen, G. Huan, and Y. Ma, Computational Methods for Multiphase Flows in Porous Media, Computational Science & Engineering, SIAM, Philadelphia, Pa, USA, 2006.
 K. Ghorayeb and A. Firoozabadi, “Numerical study of natural convection and diffusion in fractured porous media,” SPE Journal, vol. 5, no. 1, pp. 12–20, 2000. View at: Google Scholar
 H. Hoteit and A. Firoozabadi, “Multicomponent fluid flow by discontinuous Galerkin and mixed methods in unfractured and fractured media,” Water Resources Research, vol. 41, no. 11, pp. 1–15, 2005. View at: Publisher Site  Google Scholar
 H. Hoteit and A. Firoozabadi, “An efficient numerical model for incompressible twophase flow in fractured media,” Advances in Water Resources, vol. 31, no. 6, pp. 891–905, 2008. View at: Publisher Site  Google Scholar
 H. Kazemi, “Pressure transient analysis of naturally fractured reservoirs with uniform fracture distribution,” Society of Petroleum Engineers Journal, vol. 9, no. 4, pp. 451–462, 1969. View at: Google Scholar
 H. Kazeml, J. R. Gilman, and A. M. Elsharkawy, “Analytical and numerical solution of oil recovery from fractured reservoirs with empirical transfer functions,” SPE Reservoir Engineering (Society of Petroleum Engineers), vol. 7, no. 2, pp. 219–227, 1992. View at: Google Scholar
 H. Kazemi and L. S. Merrill, “Numerical simulation of water imbibition in fractured cores,” Old SPE Journal, vol. 19, no. 3, pp. 175–182, 1979. View at: Google Scholar
 S. H. Lee, C. L. Jensen, and M. F. Lough, “Efficient finitedifference model for flow in a reservoir with multiple lengthscale fractures,” SPE Journal, vol. 5, no. 3, pp. 268–275, 2000. View at: Google Scholar
 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: Google Scholar
 K. Pruess and T. N. Narasimhan, “A practical method for modeling fluid and heat flow in fractured porous media,” Society of Petroleum Engineers journal, vol. 25, no. 1, pp. 14–26, 1985. View at: Google Scholar
 S. Sarkar, M. N. Toksoz, and D. R. Burns, “Fluid flow simulation in fractured reservoirs,” in Report of the Annual Consortium Meeting, 2002. View at: Google Scholar
 L. K. Thomas, T. N. Dixon, and R. G. Pierson, “Fractured reservoir simulation,” Society of Petroleum Engineers Journal, vol. 23, no. 1, pp. 42–54, 1983. View at: Google Scholar
 J. E. Warren and P.J. Root, “The behavior of naturally fractured reservoirs,” Old SPE Journal, vol. 3, no. 3, pp. 245–255, 1963. View at: Google Scholar
 Z. Chen, G. Huan, and B. Li, “An improved IMPES method for twophase flow in porous media,” Transport in Porous Media, vol. 54, no. 3, pp. 361–376, 2004. View at: Publisher Site  Google Scholar
 K. H. Coats, “Reservoir simulation: stateoftheart,” Society of Petroleum Engineers. American Institute of Mining, Metall. and Petroleum Engineers Papers, 1982. View at: Google Scholar
 K. H. Coats, “Note on Impes and some Impesbased simulation models,” in Proceedings of the 15th SPE Reservoir Simulation Symposium, pp. 21–39, February 1999. View at: Google Scholar
 K. H. Coats, “IMPES stability: the CFL limit,” in Proceedings of the SPE Reservoir Simulation Symposium, 2001. View at: Google Scholar
 K. H. Coats, “IMPES stability: selection of stable timesteps,” SPE Journal, vol. 8, no. 2, pp. 181–187, 2003. View at: Google Scholar
 R. G. Fagin and S. CH, “A new approach to the twodimensional multiphase reservoir simulator,” Old SPE Journal, vol. 6, no. 2, pp. 175–182, 1966. View at: Google Scholar
 J. W. Sheldon, B. Zondek, and W. T. Cardwell, “Onedimensional, incompressible, noncapillary, twophase fluid flow in a porous medium,” Transactions of the American Institute of Mining and Metallurgical Engineers, vol. 216, pp. 290–296, 1959. View at: Google Scholar
 L. C. Young and R. E. Stephenson, “A generalized compositional approach for reservoir simulation,” Old SPE Journal, vol. 23, no. 5, pp. 727–742, 1983. View at: Google Scholar
 J. Kou and S. Sun, “A new treatment of capillarity to improve the stability of IMPES twophase flow formulation,” Computers and Fluids, vol. 39, no. 10, pp. 1923–1931, 2010. View at: Publisher Site  Google Scholar
 J. Kou and S. Sun, “On iterative IMPES formulation for two phase flow with capillarity in heterogeneous porous media,” International Journal of Numerical Analysis and Modeling. Series B, vol. 1, no. 1, pp. 20–40, 2010. View at: Google Scholar
 Q. Lu, A parallel multiblock/multiphysics approach for multiphase flow in porous media, Ph.D. thesis, The University of Texas at Austin, 2000.
 C. Dawson, S. Sun, and M. F. Wheeler, “Compatible algorithms for coupled flow and transport,” Computer Methods in Applied Mechanics and Engineering, vol. 193, no. 2326, pp. 2565–2580, 2004. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 V. J. Ervin, E. W. Jenkins, and S. Sun, “Coupled generalized nonlinear Stokes flow with flow through a porous medium,” SIAM Journal on Numerical Analysis, vol. 47, no. 2, pp. 929–952, 2009. View at: Publisher Site  Google Scholar
 S. Sun and J. Liu, “A locally conservative finite element method based on piecewise constant enrichment of the continuous Galerkin method,” SIAM Journal on Scientific Computing, vol. 31, no. 4, pp. 2528–2548, 2009. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. Sun and M. F. Wheeler, “A posteriori error estimation and dynamic adaptivity for symmetric discontinuous Galerkin approximations of reactive transport problems,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 78, pp. 632–652, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. Sun and M. F. Wheeler, “Anisotropic and dynamic mesh adaptation for discontinuous Galerkin methods applied to reactive transport,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 2528, pp. 3382–3405, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. Sun and M. F. Wheeler, “Projections of velocity data for the compatibility with transport,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 78, pp. 653–673, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. Sun and M. F. Wheeler, “Symmetric and nonsymmetric discontinuous Galerkin methods for reactive transport in porous media,” SIAM Journal on Numerical Analysis, vol. 43, no. 1, pp. 195–219, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 Y. Zhang, H. Wang, and T. Tang, “Simulating twophase viscoelastic flows using moving finite element methods,” Communications in Computational Physics, vol. 7, no. 2, pp. 333–349, 2010. View at: Google Scholar
 T. Belytschko and Y. Y. Lu, “Convergence and stability analyses of multitime step algorithm for parabolic systems,” Computer Methods in Applied Mechanics and Engineering, vol. 102, no. 2, pp. 179–198, 1993. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A. Gravouil and A. Combescure, “Multitimestep explicit—implicit method for nonlinear structural dynamics,” International Journal for Numerical Methods in Engineering, vol. 50, no. 1, pp. 199–225, 2000. View at: Publisher Site  Google Scholar
 M. Klisinski, “Inconsistency errors of constant velocity multitime step integration algorithms,” Computer Assisted Mechanics and Engineering Sciences, vol. 8, no. 1, pp. 121–139, 2001. View at: Google Scholar
 S. M. Bhallamudi, S. Panday, and P. S. Huyakorn, “Subtiming in fluid flow and transport simulations,” Advances in Water Resources, vol. 26, no. 5, pp. 477–489, 2003. View at: Publisher Site  Google Scholar
 Y. J. Park, E. A. Sudicky, S. Panday, J. F. Sykes, and V. Guvanasen, “Application of implicit subtime stepping to simulate flow and transport in fractured porous media,” Advances in Water Resources, vol. 31, no. 7, pp. 995–1003, 2008. View at: Publisher Site  Google Scholar
 V. Singh and S. M. Bhallamudi, “Complete hydrodynamic borderstrip irrigation model,” Journal of Irrigation and Drainage Engineering, vol. 122, no. 4, pp. 189–197, 1996. View at: Google Scholar
 V. Singh and S. M. Bhallamudi, “Hydrodynamic modeling of basin irrigation,” Journal of Irrigation and Drainage Engineering, vol. 123, no. 6, pp. 407–414, 1997. View at: Google Scholar
 P. Smolinski, T. Belytschko, and M. Neal, “Multitimestep integration using nodal partitioning,” International Journal for Numerical Methods in Engineering, vol. 26, no. 2, pp. 349–359, 1988. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 P. Smolinski, S. Sleith, and T. Belytschko, “Stability of an explicit multitime step integration algorithm for linear structural dynamics equations,” Computational Mechanics. Solids, Fluids, Engineered Materials, Aging Infrastructure, Molecular Dynamics, Heat Transfer, Manufacturing Processes, Optimization, Fracture & Integrity, vol. 18, no. 3, pp. 236–243, 1996. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 S. Sun and J. Geiser, “Multiscale discontinuous Galerkin and operatorsplitting methods for modeling subsurface flow and transport,” International Journal for Multiscale Computational Engineering, vol. 6, no. 1, pp. 87–101, 2008. View at: Publisher Site  Google Scholar
 J. E. VanderKwaak, Numerical simulation of flow and chemical transport in integrated surfacesubsurface hydrologic systems, Ph.D. thesis, University of Waterloo, 1999.
 P. A. Forsyth, Jr. and P. H. Sammon, “Quadratic convergence for cellcentered grids,” Applied Numerical Mathematics, vol. 4, no. 5, pp. 377–394, 1988. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. E. P. Monteagudo and A. Firoozabadi, “Controlvolume method for numerical simulation of twophase immiscible flow in two and threedimensional discretefractured media,” Water Resources Research, vol. 40, no. 7, pp. W074051–W0740520, 2004. View at: Google Scholar
Copyright
Copyright © 2011 Jisheng Kou 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.