Abstract

The temporal discretization scheme is one important ingredient of efficient simulator for two-phase flow in the fractured porous media. The application of single-scale 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 multi-scale time splitting strategy to simulate multi-scale multi-physics processes of two-phase flow in fractured porous media. We use the multi-scale time schemes for both the pressure and saturation equations; that is, a large time-step size is employed for the matrix domain, along with a small time-step 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 two-phase 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 matrix-fracture 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 single-porosity model, the dual-porosity/dual-permeability model, and the discrete-fracture model. In the single-porosity 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 dual-porosity/dual-permeability 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 single-porosity model. The discrete fracture model is shown to be is numerically superior to the single-porosity model and the dual-porosity models [6, 7]. This model will be used in this paper.

The implicit pressure explicit saturation (IMPES) method [17–24] is a popular time-stepping 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 time-step 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 two-phase 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 time-step 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 time-step size to match the rapid changes in the fracture network if using the single-scale time-stepping. To resolve this problem, we should use the multiscale time-stepping 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 fine-grid locations. We fist carry out the computations for all coarse grid-blocks by using a time-step size that satisfies the coarse grid-block stability constraint. Once the solutions in a coarsely discretized region are obtained, the unknowns within the fine grids are solved by using a smaller time-step size that matches the fine grid stability constraints. In the subtiming methodology for implicit-type time-stepping schemes proposed in [39], the fully implicit schemes with multiple time-step sizes are utilized to the different portions of a domain, which have different accuracy requirements. This technique has been applied to the time-dependent flow and transport simulations in fractured porous media [40], which uses the implicit subtime-stepping in the locations within and near the fractures. It is demonstrated that the implicit multiple time-stepping is capable to significantly improve the accuracy and efficiency of numerical simulation. The hybrid implicit-explicit approach [46] treats implicitly the portions of the domain with the particular stability requirements, along with the explicit time-stepping 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 time-step.

In this paper, we present a multiscale time splitting strategy for two-phase 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 two-phase 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 time-step size in the matrix domain, along with a small time-step size being applied in the fracture system. The time-step for the pressure in the fractures is further partitioned to a few subtime-steps, in which the saturation will be updated explicitly. The updating process of the saturation employs two-scale time-stepping; a coarse time-stepping 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 two-phase incompressible flow model, along with the discrete-fracture model. In Section 3, we describe our multiscale time splitting strategy. Here, the cell-centered 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 Two-Phase Incompressible Flow

2.1. Governing Equations of Two-Phase Flow

We now describe two-phase 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𝜙𝜕𝑆𝛼𝜕𝑡+∇⋅𝐮𝛼=𝑞𝛼,𝛼=𝑤,𝑛,(2.1) 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𝐮𝛼𝑘=−𝑟𝛼𝜇𝛼𝐊∇𝑝𝛼+𝜌𝛼𝑔∇𝑧,𝛼=𝑤,𝑛,(2.2) 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𝑆𝑤+𝑆𝑛=1.(2.3) The difference between the nonwetting phase and wetting phase pressures is described by the capillary pressure𝑝𝑐𝑆𝑤=𝑝𝑛−𝑝𝑤.(2.4)

In this paper, we use the two-phase flow formulation [1, 7, 48] given by𝐮∇⋅𝑎+𝐮𝑐≡−∇⋅𝜆𝑡𝐊∇Φ𝑤−∇⋅𝜆𝑛𝐊∇Φ𝑐=𝑞𝑡,(2.5)𝜙𝜕𝑆𝑤𝜕𝑡−𝑞𝑤𝑓=−∇⋅𝑤𝐮𝑎≡∇⋅𝜆𝑤𝐊∇Φ𝑤,(2.6) where 𝑞𝑡=𝑞𝑤+𝑞𝑛, 𝜆𝛼=𝑘𝑟𝛼/𝜇𝛼, 𝜆𝑡=𝜆𝑤+𝜆𝑛, 𝑓𝑤=𝜆𝑤/𝜆𝑡, Φ𝛼=𝑝𝛼+𝜌𝑔𝑧, Φ𝑐=𝑝𝑐+(𝜌𝑛−𝜌𝑤)𝑔𝑧, 𝐮𝑎=−𝜆𝑡𝐊∇Φ𝑤, and 𝐮𝑐=−𝜆𝑛𝐊∇Φ𝑐. From the above definitions, the wetting-phase velocity is expressed as𝐮𝑤=𝑓𝑤𝐮𝑎.(2.7)

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𝑝𝑤or𝑝𝑛=𝑝𝐷onΓ𝐷,𝐮𝑎+𝐮𝑐⋅𝐧=𝐮𝑁onΓ𝑁,(2.8) 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𝑆𝑤or𝑆𝑛=𝑆𝑁onΓ𝑁.(2.9) The initial saturation of the wetting phase is given by𝑆𝑤=𝑆0𝑤inΩ.(2.10)

2.2. Discrete-Fracture Model

Here, we employ the discrete-fracture 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 (𝑛−1)-dimensional. This treatment is capable to considerably improve the computational efficiency since it can remove the length-scale contrast resulting from the explicit representation of the fracture aperture as in the single-porosity model. The entire domain is decomposed into two parts: the matrix Ω𝑚 (𝑛-dimensional) and fracture Ω𝑓 ((𝑛−1)-dimensional), and the fractures are surrounded by the matrix blocks. The pressure equation in the matrix domain is given by−∇⋅𝜆𝑡,𝑚𝐊𝑚∇Φ𝑤,𝑚−∇⋅𝜆𝑛,𝑚𝐊𝑚∇Φ𝑐,𝑚=𝑞𝑡,𝑚,(2.11) where the subscript 𝑚 represents the matrix domain. Equation (2.11) is subject to the matrix-fracture interface conditionΦ𝑤,𝑚=Φ𝑤,𝑓,Φ𝑐,𝑚=Φ𝑐,𝑓,on𝜕Ω𝑚∩Ω𝑓.(2.12) 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−∇⋅𝜆𝑡,𝑓𝐊𝑓∇Φ𝑤,𝑓−∇⋅𝜆𝑛,𝑓𝐊𝑓∇Φ𝑐,𝑓=𝑞𝑡,𝑓+𝑄𝑡,𝑓,(2.13) where 𝑄𝑡,𝑓 is the mass transfer across the matrix-fracture interfaces. Similarly, we can express the saturation equation in the matrix domain as𝜙𝑚𝜕𝑆𝑤,𝑚𝜕𝑡−∇⋅𝜆𝑤,𝑚𝐊𝑚∇Φ𝑤,𝑚=𝑞𝑤,𝑚,(2.14) with the interface condition𝑆𝑤,𝑚=𝑆𝑤,𝑓,on𝜕Ω𝑚∩Ω𝑓.(2.15) In accordance, the saturation equation in the fracture system is given by𝜙𝑓𝜕𝑆𝑤,𝑓𝜕𝑡−∇⋅𝜆𝑤,𝑓𝐊𝑓∇Φ𝑤,𝑓=𝑞𝑤,𝑓+𝑄𝑤,𝑓,(2.16) where 𝑄𝑤,𝑓 represents the mass transfer across the matrix-fracture 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𝜙𝑚𝜕𝑆𝑤,𝑚𝑓𝜕𝑡+∇⋅𝑤,𝑚𝐮𝑎,𝑚=𝑞𝑤,𝑚,(2.17) and in the fractures𝜙𝑓𝜕𝑆𝑤,𝑓𝑓𝜕𝑡+∇⋅𝑤,𝑓𝐮𝑎,𝑓=𝑞𝑤,𝑓+𝑄𝑤,𝑓,(2.18) where 𝑄𝑤,𝑓 is the same to 𝑄𝑤,𝑓, but with different expressions in the spatial discretization.

3. Multi-Scale Time Splitting Strategy

In this paper, we use the cell-centered 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 two-dimensional space. The extension to three-dimensional case is straightforward. We assume that the porous media are isotropic and the absolute permeability is a scalar function.

3.1. Multi-Scale Time Splitting Approach for the Pressure Equation with Implicit Capillary Pressure

We firstly divide the total time interval [0,𝑇] into 𝑁𝑝,𝑚 time-steps as 0=𝑡0<𝑡1<⋯<𝑡𝑁𝑝,𝑚=𝑇. This time division is used for the pressure in the matrix domain. Define the time-step length Δ𝑡𝑖=𝑡𝑖+1−𝑡𝑖. Since the pressure in the fractures varies more rapidly than the matrix domain, we use a smaller time-step size for the pressure in the fractures; that is, we partition each subinterval (𝑡𝑖,𝑡𝑖+1] into 𝑁𝑝,𝑓 subsubintervals as (𝑡𝑖,𝑡𝑖+1⋃]=𝑁𝑝,𝑓−1𝑗=0(𝑡𝑖,𝑗,𝑡𝑖,𝑗+1], where 𝑡𝑖,0=𝑡𝑖 and 𝑡𝑖,𝑁𝑝,𝑓=𝑡𝑖+1, and denote the subsubinterval length by Δ𝑡𝑖,𝑗=𝑡𝑖,𝑗+1−𝑡𝑖,𝑗. 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 one-dimensional 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 time-step 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−∇⋅𝜆𝑡𝑆𝑖𝑤,𝑚𝐊𝑚∇Φ𝑖+1𝑤,𝑚−∇⋅𝜆𝑛𝑆𝑖𝑤,𝑚𝐊𝑚∇(1−𝜔)Φ𝑖𝑐,𝑚+𝜔Φ𝑖+1𝑐,𝑚=𝑞𝑖+1𝑡,𝑚,(3.1) where the superscript 𝑖 represents the time-step and 𝜔 is the temporal discretization parameter. Apparently, it becomes the form of IMPES if 𝜔=0. Here, we always take 𝜔≥1/2; that is, the capillarity is implicitly treated. It is similar to obtain the form in the fracture (referred to by the subscript 𝑓) as−∇⋅𝜆𝑡𝑆𝑖𝑤,𝑓𝐊𝑓∇Φ𝑖+1𝑤,𝑓−∇⋅𝜆𝑛𝑆𝑖𝑤,𝑓𝐊𝑓∇(1−𝜔)Φ𝑖𝑐,𝑓+𝜔Φ𝑖+1𝑐,𝑓=𝑞𝑖+1𝑡,𝑓+𝑄𝑖+1𝑡,𝑓.(3.2)

Let 𝑒 indicate the boundary of the matrix gridcells. We use the CCFD scheme applied to (3.1) and obtain𝑒∈𝜕𝐾𝐹𝑖+1𝑎,𝑚,𝑒+𝑒∈𝜕𝐾(1−𝜔)𝐹𝑖𝑐,𝑚,𝑒+𝜔𝐹𝑖+1𝑐,𝑚,𝑒=𝑞𝑖+1𝑡,𝑚,𝐾||𝐾||,(3.3) 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 𝐹𝑖+1𝑎,𝑚,𝑒=−|𝑒|𝛽𝑖𝑡,𝑒Φ𝑖+1𝑤,𝑚,𝐾′−Φ𝑖+1𝑤,𝑚,𝐾𝑑𝐾,𝐾′,𝐹𝑖𝑐,𝑚,𝑒=−|𝑒|𝛽𝑖𝑛,𝑒Φ𝑖𝑐,𝑚,𝐾′−Φ𝑖𝑐,𝑚,𝐾𝑑𝐾,𝐾′,𝐹𝑖+1𝑐,𝑚,𝑒=−|𝑒|𝛽𝑖𝑛,𝑒Φ𝑖+1𝑐,𝑚,𝐾′−Φ𝑖+1𝑐,𝑚,𝐾𝑑𝐾,𝐾′,(3.4) 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𝛽𝑖𝑡,𝑒=𝑑𝐾,𝐾′𝜆𝑖𝑡,𝑚,𝐾𝜆𝑖𝑡,𝑚,𝐾′𝐊𝑚,𝐾𝐊𝑚,𝐾′𝑑𝐾,𝑒𝜆𝑖𝑡,𝑚,𝐾′𝐊𝑚,𝐾′+𝑑𝐾′,𝑒𝜆𝑖𝑡,𝑚,𝐾𝐊𝑚,𝐾,𝛽𝑖𝑛,𝑒=𝑑𝐾,𝐾′𝜆𝑖𝑛,𝑚,𝐾𝜆𝑖𝑛,𝑚,𝐾′𝐊𝑚,𝐾𝐊𝑚,𝐾′𝑑𝐾,𝑒𝜆𝑖𝑛,𝑚,𝐾′𝐊𝑚,𝐾′+𝑑𝐾′,𝑒𝜆𝑖𝑛,𝑚,𝐾𝐊𝑚,𝐾,(3.5) 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𝐹𝑖+1𝑎,𝑚,𝑒≡𝐹𝑖+1𝑎,𝑚,𝑒,𝐾=−|𝑒|𝛽𝑖𝑡,𝑚𝑓,𝑒Φ𝑖+1𝑤,𝑓,𝑒−Φ𝑖+1𝑤,𝑚,𝐾𝑑𝐾,𝑒,+𝜀/2(3.6)𝐹𝑖𝑐,𝑚,𝑒≡𝐹𝑖𝑐,𝑚,𝑒,𝐾=−|𝑒|𝛽𝑖𝑛,𝑚𝑓,𝑒Φ𝑖𝑐,𝑓,𝑒−Φ𝑖𝑐,𝑚,𝐾𝑑𝐾,𝑒,+𝜀/2(3.7)𝐹𝑖+1𝑐,𝑚,𝑒≡𝐹𝑖+1𝑐,𝑚,𝑒,𝐾=−|𝑒|𝛽𝑖𝑛,𝑚𝑓,𝑒Φ𝑖+1𝑐,𝑓,𝑒−Φ𝑖+1𝑐,𝑚,𝐾𝑑𝐾,𝑒,+𝜀/2(3.8) where 𝛽𝑖𝑡,𝑚𝑓,𝑒 and 𝛽𝑖𝑛,𝑚𝑓,𝑒 are given by𝛽𝑖𝑡,𝑚𝑓,𝑒=𝑑𝐾,𝑒𝜆+𝜀/2𝑖𝑡,𝑓,𝐾𝜆𝑖𝑡,𝑚,𝐾𝐊𝑓,𝐾𝐊𝑚,𝐾(𝜀/2)𝜆𝑖𝑡,𝑚,𝐾𝐊𝑚,𝐾+𝑑𝐾,𝑒𝜆𝑖𝑡,𝑓,𝐾𝐊𝑓,𝐾,𝛽𝑖𝑛,𝑚𝑓,𝑒=𝑑𝐾,𝑒𝜆+𝜀/2𝑖𝑛,𝑓,𝐾𝜆𝑖𝑛,𝑚,𝐾𝐊𝑓,𝐾𝐊𝑚,𝐾(𝜀/2)𝜆𝑖𝑛,𝑚,𝐾𝐊𝑚,𝐾+𝑑𝐾,𝑒𝜆𝑖𝑛,𝑓,𝐾𝐊𝑓,𝐾.(3.9) 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 reduced-dimensional fracture system. Let 𝑒 be a gridcell of the fracture network. The CCFD scheme is applied to (3.2), and then we have𝛾∈𝜕𝑒𝐹𝑖+1𝑎,𝑓,𝛾+𝛾∈𝜕𝑒(1−𝜔)𝐹𝑖𝑐,𝑓,𝛾+𝜔𝐹𝑖+1𝑐,𝑓,𝛾=𝑞𝑖+1𝑡,𝑓,𝑒|𝑒|+𝑄𝑖+1𝑡,𝑓,𝑒|𝑒|,(3.10) 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 fracture-fracture fluxes in (3.10), but the matrix-fracture 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 matrix-fracture interface 𝑒 is given by𝑄𝑖+1𝑡,𝑓,𝑒𝑄=−𝑖+1𝑡,𝑓,𝑒,𝐾+𝑄𝑖+1𝑡,𝑓,𝑒,𝐾′𝜀,𝑄𝑖+1𝑡,𝑓,𝑒,𝐾=𝐹𝑖+1𝑎,𝑚,𝑒,𝐾+(1−𝜔)𝐹𝑖𝑐,𝑚,𝑒,𝐾+𝜔𝐹𝑖+1𝑐,𝑚,𝑒,𝐾,𝑄𝑖+1𝑡,𝑓,𝑒,𝐾′=𝐹𝑖+1𝑎,𝑚,𝑒,𝐾′+(1−𝜔)𝐹𝑖𝑐,𝑚,𝑒,𝐾′+𝜔𝐹𝑖+1𝑐,𝑚,𝑒,𝐾′,(3.11) where 𝐹𝑖+1𝑎,𝑚,𝑒, 𝐹𝑖𝑐,𝑚,𝑒, and 𝐹𝑖+1𝑐,𝑚,𝑒 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𝑒∈Λ𝛾𝐹𝑖+1𝑎,𝑓,𝛾,𝑒+(1−𝜔)𝐹𝑖𝑐,𝑓,𝛾,𝑒+𝜔𝐹𝑖+1𝑐,𝑓,𝛾,𝑒=0,(3.12) where 𝐹𝑎,𝑓,𝛾,𝑒=𝐹𝑎,𝑓,𝛾|𝛾∈𝑒 and 𝐹𝑐,𝑓,𝛾,𝑒=𝐹𝑐,𝑓,𝛾|𝛾∈𝑒.

Combining the formulations in the matrix domain and fracture network, we obtain the discretization of total mass conservation equation given by𝐀𝑖𝑎,𝑚,𝑚𝐀𝑖𝑎,𝑚,𝑓𝐀𝑖𝑎,𝑓,𝑚𝐀𝑖𝑎,𝑓,𝑓Φ𝑖+1𝑤,𝑚Φ𝑖+1𝑤,𝑓+𝐀𝑖𝑐,𝑚,𝑚𝐀𝑖𝑐,𝑚,𝑓𝐀𝑖𝑐,𝑓,𝑚𝐀𝑖𝑐,𝑓,𝑓×Φ(1−𝜔)𝑖𝑐,𝑚Φ𝑖𝑐,𝑓Φ+𝜔𝑖+1𝑐,𝑚Φ𝑖+1𝑐,𝑓=𝐐𝑖+1𝑎𝑐,𝑚𝐐𝑖+1𝑎𝑐,𝑓.(3.13) The terms 𝐀𝑖𝑎,𝑚,𝑓, 𝐀𝑖𝑎,𝑓,𝑚, 𝐀𝑖𝑐,𝑚,𝑓, and 𝐀𝑖𝑐,𝑓,𝑚 of (3.13) arise from the interconnections of the matrix blocks and fracture system. Since 𝜔>0 and 𝐒𝑤𝑖+1 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 Φ𝑐(𝐒𝑤𝑖+1) at 𝐒𝑖𝑤 byΦ𝑐𝐒𝑤𝑖+1≃Φ𝑐𝐒𝑖𝑤+Φ𝑐𝐒𝑖𝑤𝐒𝑤𝑖+1−𝐒𝑖𝑤.(3.14) Furthermore, we employ the backward Euler time discretization for the saturation equation both in the matrix domain𝜙𝑚𝑆𝑖+1𝑤,𝑚−𝑆𝑖𝑤,𝑚Δ𝑡𝑖−∇⋅𝜆𝑤𝑆𝑖𝑤,𝑚𝐊𝑚∇Φ𝑖+1𝑤,𝑚=𝑞𝑖+1𝑤,𝑚,(3.15) and in the fracture network𝜙𝑓𝑆𝑖+1𝑤,𝑓−𝑆𝑖𝑤,𝑓Δ𝑡𝑖−∇⋅𝜆𝑤𝑆𝑖𝑤,𝑓𝐊𝑓∇Φ𝑖+1𝑤,𝑓=𝑞𝑖+1𝑤,𝑓+𝑄𝑖+1𝑤,𝑓.(3.16) It is analogous to the derivation of the pressure equation (3.13), and we approximate (3.15) and (3.16) by the CCFD method as1Δ𝑡𝑖𝐌𝑚𝐌𝑓𝐒𝑖+1𝑤,𝑚−𝐒𝑖𝑤,𝑚𝐒𝑖+1𝑤,𝑓−𝐒𝑖𝑤,𝑓+𝐀𝑖𝑤,𝑚,𝑚𝐀𝑖𝑤,𝑚,𝑓𝐀𝑖𝑤,𝑓,𝑚𝐀𝑖𝑤,𝑓,𝑓Φ𝑖+1𝑤,𝑚Φ𝑖+1𝑤,𝑓=𝐐𝑖+1𝑤,𝑚𝐐𝑖+1𝑤,𝑓.(3.17) Substituting (3.14) and (3.17) into (3.13), we obtain the coupled pressure equation𝐀𝑖𝑡Φ𝑖+1𝑤,𝑚Φ𝑖+1𝑤,𝑓=𝐐𝑖𝑡,(3.18) where𝐀𝑖𝑡=𝐀𝑖𝑎,𝑚,𝑚𝐀𝑖𝑎,𝑚,𝑓𝐀𝑖𝑎,𝑓,𝑚𝐀𝑖𝑎,𝑓,𝑓−𝜔Δ𝑡𝑖𝐀𝑖𝑐,𝑚,𝑚𝐀𝑖𝑐,𝑚,𝑓𝐀𝑖𝑐,𝑓,𝑚𝐀𝑖𝑐,𝑓,𝑓⎡⎢⎢⎣Φ𝑐𝐒𝑖𝑤,𝑚Φ𝑐𝐒𝑖𝑤,𝑓⎤⎥⎥⎦×𝐌𝑚−1𝐌𝑓−1𝐀𝑖𝑤,𝑚,𝑚𝐀𝑖𝑤,𝑚,𝑓𝐀𝑖𝑤,𝑓,𝑚𝐀𝑖𝑤,𝑓,𝑓,𝐐𝑖𝑡=𝐐𝑖+1𝑎𝑐,𝑚𝐐𝑖+1𝑎𝑐,𝑓−𝐀𝑖𝑐,𝑚,𝑚𝐀𝑖𝑐,𝑚,𝑓𝐀𝑖𝑐,𝑓,𝑚𝐀𝑖𝑐,𝑓,𝑓⎡⎢⎢⎣Φ𝑐𝐒𝑖𝑤,𝑚Φ𝑐𝐒𝑖𝑤,𝑓⎤⎥⎥⎦−𝜔Δ𝑡𝑖𝐀𝑖𝑐,𝑚,𝑚𝐀𝑖𝑐,𝑚,𝑓𝐀𝑖𝑐,𝑓,𝑚𝐀𝑖𝑐,𝑓,𝑓×⎡⎢⎢⎣Φ𝑐𝐒𝑖𝑤,𝑚Φ𝑐𝐒𝑖𝑤,𝑓⎤⎥⎥⎦𝐌𝑚−1𝐌𝑓−1𝐐𝑖+1𝑤,𝑚𝐐𝑖+1𝑤,𝑓.(3.19)

Equation (3.18), along with (3.19), is utilized to compute the matrix pressure in a large time-step. In order to match the rapid changes of the pressure in the fracture network, we use a small time-step size in the fracture system. For each time substep, we consider the pressure equation in the fractures as−∇⋅𝜆𝑡𝑆𝑖,𝑗𝑤,𝑓𝐊𝑓∇Φ𝑖,𝑗+1𝑤,𝑓−∇⋅𝜆𝑛𝑆𝑖,𝑗𝑤,𝑓𝐊𝑓∇(1−𝜔)Φ𝑖,𝑗𝑐,𝑓+𝜔Φ𝑖,𝑗+1𝑐,𝑓=𝑞𝑖,𝑗+1𝑡,𝑓+𝑄𝑖,𝑗+1𝑡,𝑓.(3.20) It is obtained by using the CCFD scheme to (3.20) that𝛾∈𝜕𝑒𝐹𝑖,𝑗+1𝑎,𝑓,𝛾+𝛾∈𝜕𝑒(1−𝜔)𝐹𝑖,𝑗𝑐,𝑓,𝛾+𝜔𝐹𝑖,𝑗+1𝑐,𝑓,𝛾=𝑞𝑖,𝑗+1𝑡,𝑓,𝑒|𝑒|+𝑄𝑖,𝑗+1𝑡,𝑓,𝑒|𝑒|,(3.21) 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 matrix-fracture transfer is given by𝑄𝑖,𝑗+1𝑡,𝑓,𝑒𝑄=−𝑖,𝑗+1𝑡,𝑓,𝑒,𝐾+𝑄𝑖,𝑗+1𝑡,𝑓,𝑒,𝐾′𝜀,𝑄𝑖,𝑗+1𝑡,𝑓,𝑒,𝐾=𝐹𝑖,𝑗+1𝑎,𝑚,𝑒,𝐾+(1−𝜔)𝐹𝑖,𝑗𝑐,𝑚,𝑒,𝐾+𝜔𝐹𝑖,𝑗+1𝑐,𝑚,𝑒,𝐾,𝑄𝑖,𝑗+1𝑡,𝑓,𝑒,𝐾′=𝐹𝑖,𝑗+1𝑎,𝑚,𝑒,𝐾′+(1−𝜔)𝐹𝑖,𝑗𝑐,𝑚,𝑒,𝐾′+𝜔𝐹𝑖,𝑗+1𝑐,𝑚,𝑒,𝐾′,(3.22) where𝐹𝑖,𝑗+1𝑎,𝑚,𝑒,𝐾≡𝐹𝑖,𝑗+1𝑎,𝑚,𝑒=−|𝑒|𝛽𝑖,𝑗𝑡,𝑚𝑓,𝑒Φ𝑖,𝑗+1𝑤,𝑓,𝑒−Φ𝑖+1𝑤,𝑚,𝐾𝑑𝐾,𝑒,+𝜀/2(3.23)𝐹𝑖,𝑗𝑐,𝑚,𝑒,𝐾≡𝐹𝑖,𝑗𝑐,𝑚,𝑒=−|𝑒|𝛽𝑖,𝑗𝑐,𝑚𝑓,𝑒Φ𝑖,𝑗𝑐,𝑓,𝑒−Φ𝑖,𝑗𝑐,𝑚,𝐾𝑑𝐾,𝑒,𝐹+𝜀/2𝑖,𝑗+1𝑐,𝑚,𝑒,𝐾≡𝐹𝑖,𝑗+1𝑐,𝑚,𝑒=−|𝑒|𝛽𝑖,𝑗𝑐,𝑚𝑓,𝑒Φ𝑖,𝑗+1𝑐,𝑓,𝑒−Φ𝑖,𝑗+1𝑐,𝑚,𝐾𝑑𝐾,𝑒,+𝜀/2(3.24) along with𝛽𝑖,𝑗𝑡,𝑚𝑓,𝑒=𝑑𝐾,𝑒𝜆+𝜀/2𝑖,𝑗𝑡,𝑓,𝐾𝜆𝑖,𝑗𝑡,𝑚,𝐾𝐊𝑓,𝐾𝐊𝑚,𝐾(𝜀/2)𝜆𝑖,𝑗𝑡,𝑚,𝐾𝐊𝑚,𝐾+𝑑𝐾,𝑒𝜆𝑖,𝑗𝑡,𝑓,𝐾𝐊𝑓,𝐾,𝛽𝑖,𝑗𝑐,𝑚𝑓,𝑒=𝑑𝐾,𝑒𝜆+𝜀/2𝑖,𝑗𝑐,𝑓,𝐾𝜆𝑖,𝑗𝑐,𝑚,𝐾𝐊𝑓,𝐾𝐊𝑚,𝐾(𝜀/2)𝜆𝑖,𝑗𝑐,𝑚,𝐾𝐊𝑚,𝐾+𝑑𝐾,𝑒𝜆𝑖,𝑗𝑐,𝑓,𝐾𝐊𝑓,𝐾.(3.25) The treatment of the interface connecting multiple fractures is similar to (3.12). Therefore, in the time-step used for the pressure in the fractures, we solve the following pressure equation𝐀𝑖,𝑗𝑡,𝑓Φ𝑖,𝑗+1𝑤,𝑓=𝐐𝑖,𝑗𝑡,𝑓,(3.26) where 𝐀𝑖,𝑗𝑡,𝑓=𝐀𝑖,𝑗𝑎,𝑓,𝑓−𝜔Δ𝑡𝑖,𝑗𝐀𝑖,𝑗𝑐,𝑓,𝑓Φ𝑐𝐒𝑖,𝑗𝑤,𝑓𝐌𝑓−1𝐀𝑖,𝑗𝑤,𝑓,𝑓+𝐀𝑖,𝑗𝑐,𝑓,𝑚Φ𝑐𝐒𝑖,𝑗𝑤,𝑚𝐌𝑚−1𝐀𝑖,𝑗𝑤,𝑚,𝑓,𝐐𝑖,𝑗𝑡,𝑓=𝐐𝑖,𝑗+1𝑎𝑐,𝑓−𝐀𝑖,𝑗𝑎,𝑓,𝑚Φ𝑖+1𝑤,𝑚+𝜔Δ𝑡𝑖,𝑗𝐀𝑖,𝑗𝑐,𝑓,𝑚Φ𝑐𝐒𝑖,𝑗𝑤,𝑚𝐌𝑚−1𝐀𝑖,𝑗𝑤,𝑚,𝑚+𝐀𝑖,𝑗𝑐,𝑓,𝑓Φ𝑐𝐒𝑖,𝑗𝑤,𝑓𝐌𝑓−1𝐀𝑖,𝑗𝑤,𝑓,𝑚Φ𝑖+1𝑤,𝑚−𝐀𝑖,𝑗𝑐,𝑓,𝑚Φ𝑐𝐒𝑖,𝑗𝑤,𝑚−𝐀𝑖,𝑗𝑐,𝑓,𝑓Φ𝑐𝐒𝑖,𝑗𝑤,𝑓−𝜔Δ𝑡𝑖,𝑗𝐀𝑖,𝑗𝑐,𝑓,𝑚Φ𝑐𝐒𝑖,𝑗𝑤,𝑚𝐌𝑚−1𝐐𝑖+1𝑤,𝑚+𝐀𝑖,𝑗𝑐,𝑓,𝑓Φ𝑐𝐒𝑖,𝑗𝑤,𝑓𝐌𝑓−1𝐐𝑖+1𝑤,𝑓.(3.27)

In the above equations, the inverse of 𝐌𝑚 and 𝐌𝑓 is not expensive as they are diagonal. In the time-step (𝑡𝑖,𝑡𝑖+1), we solve the linear system (3.18) implicitly to obtain the pressure Φ𝑖+1𝑤,𝑚, and then in the time substep (𝑡𝑖,𝑗,𝑡𝑖,𝑗+1) compute Φ𝑖,𝑗+1𝑤,𝑓 by solving (3.26).

Once the pressures Φ𝑖+1𝑤,𝑚 and Φ𝑖,𝑗+1𝑤,𝑓 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 𝑒̸⊆Ω𝑓,𝐹𝑖,𝑗+1𝑎,𝑚,𝑒=−|𝑒|𝛽𝑖,𝑗𝑡,𝑒Φ𝑖+1𝑤,𝑚,𝐾′−Φ𝑖+1𝑤,𝑚,𝐾𝑑𝐾,𝐾′,(3.28) where𝛽𝑖,𝑗𝑡,𝑒=𝑑𝐾,𝐾′𝜆𝑖,𝑗𝑡,𝑚,𝐾𝜆𝑖,𝑗𝑡,𝑚,𝐾′𝐊𝑚,𝐾𝐊𝑚,𝐾′𝑑𝐾,𝑒𝜆𝑖,𝑗𝑡,𝑚,𝐾′𝐊𝑚,𝐾′+𝑑𝐾′,𝑒𝜆𝑖,𝑗𝑡,𝑚,𝐾𝐊𝑚,𝐾.(3.29) 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. Multi-Scale Explicit Time-Stepping for the Saturation Equation

As the pressure computed by the above-described method has coupled with the saturation, the forward Euler time discretization is used for computing the saturation equation of the wetting phase in each time-step used for the pressure in the fracture. According to the concept proposed in [17, 24, 27], the time-step 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 time-step (𝑡𝑖,𝑗,𝑡𝑖,𝑗+1] for the pressure in the fracture network into 𝑁𝑠,𝑚 time-steps as (𝑡𝑖,𝑗,𝑡𝑖,𝑗+1⋃]=𝑁𝑠,𝑚−1𝑘=0(𝑡𝑖,𝑗,𝑘,𝑡𝑖,𝑗,𝑘+1], where 𝑡𝑖,𝑗,0=𝑡𝑖,𝑗 and 𝑡𝑖,𝑗,𝑁𝑠,𝑚=𝑡𝑖,𝑗+1. 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 time-step size is required for the saturation in the fracture system; that is, we partition (𝑡𝑖,𝑗,𝑘,𝑡𝑖,𝑗,𝑘+1] into 𝑁𝑠,𝑓 time substeps as (𝑡𝑖,𝑗,𝑘,𝑡𝑖,𝑗,𝑘+1⋃]=𝑁𝑠,𝑓−1𝑙=0(𝑡𝑖,𝑗,𝑘,𝑙,𝑡𝑖,𝑗,𝑘,𝑙+1], where 𝑡𝑖,𝑗,𝑘,0=𝑡𝑖,𝑗,𝑘 and 𝑡𝑖,𝑗,𝑘,𝑁𝑠,𝑓=𝑡𝑖,𝑗,𝑘+1. Denote Δ𝑡𝑖,𝑗,𝑘=𝑡𝑖,𝑗,𝑘+1−𝑡𝑖,𝑗,𝑘 and Δ𝑡𝑖,𝑗,𝑘,𝑙=𝑡𝑖,𝑗,𝑘,𝑙+1−𝑡𝑖,𝑗,𝑘,𝑙. The explicit scheme is employed for the saturation equation both in the matrix domain𝜙𝑚𝑆𝑖,𝑗,𝑘+1𝑤,𝑚−𝑆𝑖,𝑗,𝑘𝑤,𝑚Δ𝑡𝑖,𝑗,𝑘𝑓+∇⋅𝑖,𝑗,𝑘𝑤,𝑚𝐮𝑖,𝑗+1𝑎,𝑚=𝑞𝑖,𝑗,𝑘𝑤,𝑚,(3.30) and in the fracture network𝜙𝑓𝑆𝑖,𝑗,𝑘,𝑙+1𝑤,𝑓−𝑆𝑖,𝑗,𝑘,𝑙𝑤,𝑓Δ𝑡𝑖,𝑗,𝑘,𝑙𝑓+∇⋅𝑖,𝑗,𝑘,𝑙𝑤,𝑓𝐮𝑖,𝑗+1𝑎,𝑓=𝑞𝑖,𝑗,𝑘,𝑙𝑤,𝑓+𝑄𝑖,𝑗,𝑘𝑤,𝑓.(3.31) We use the upwind CCFD method to discretize the saturation equation (3.30)||𝐾||𝜙𝑚,𝐾𝑆𝑖,𝑗,𝑘+1𝑤,𝑚,𝐾−𝑆𝑖,𝑗,𝑘𝑤,𝑚,𝐾Δ𝑡𝑖,𝑗,𝑘+𝑒∈𝜕𝐾𝑓𝑖,𝑗,𝑘𝑤,𝑒𝐹𝑖,𝑗+1𝑎,𝑚,𝑒=𝑞𝑖,𝑗,𝑘𝑤,𝑚,𝐾||𝐾||.(3.32) Let 𝑒 be the interface between the matrix gridcells 𝐾 and 𝐾; that is, 𝑒=𝜕𝐾∩𝜕𝐾. If 𝑒̸⊆Ω𝑓, the term 𝑓𝑤,𝑒 in (3.32) is determined by𝑓𝑖,𝑗,𝑘𝑤,𝑒=𝑓𝑖,𝑗,𝑘𝑤,𝑚,𝐾,𝐹𝑖,𝑗+1𝑎,𝑚,𝑒𝑓>0,𝑖,𝑗,𝑘𝑤,𝑚,𝐾′,𝐹𝑖,𝑗+1𝑎,𝑚,𝑒<0.(3.33) If 𝑒⊆Ω𝑓 is a gridcell of the fracture network, the term 𝑓𝑤,𝑒 in (3.32) is determined by𝑓𝑖,𝑗,𝑘𝑤,𝑒=𝑓𝑖,𝑗,𝑘𝑤,𝑚,𝐾,𝐹𝑖,𝑗+1𝑎,𝑚,𝑒𝑓>0,𝑖,𝑗,𝑘𝑤,𝑓,𝑒,𝐹𝑖,𝑗+1𝑎,𝑚,𝑒<0,(3.34) where𝑓𝑖,𝑗,𝑘𝑤,𝑓,𝑒=1𝑁𝑁𝑠,𝑓𝑠,𝑓−1𝑙=0𝑓𝑖,𝑗,𝑘,𝑙+1𝑤,𝑓,𝑒+𝑓𝑖,𝑗,𝑘,𝑙𝑤,𝑓,𝑒2.(3.35)

It is analogous to discretize the saturation equation in the fracture system as|𝑒|𝜙𝑓,𝑒𝑆𝑖,𝑗,𝑘,𝑙+1𝑤,𝑓,𝑒−𝑆𝑖,𝑗,𝑘,𝑙𝑤,𝑓,𝑒Δ𝑡𝑖,𝑗,𝑘,𝑙+𝛾∈𝜕𝑒𝑓𝑖,𝑗,𝑘,𝑙𝑤,𝛾𝐹𝑖,𝑗+1𝑎,𝑓,𝛾=𝑞𝑖,𝑗,𝑘,𝑙𝑤,𝑓,𝛾|𝑒|+𝑄𝑖,𝑗,𝑘,𝑙𝑤,𝑓,𝑒|𝑒|,(3.36) where 𝑒 is a gridcell of the fracture network. Let 𝛾=𝜕𝑒∩𝜕𝑒 where 𝑒 and 𝑒 are the fracture gridcells, then we have𝑓𝑖,𝑗,𝑘,𝑙𝑤,𝛾=𝑓𝑖,𝑗,𝑘,𝑙𝑤,𝑓,𝑒,𝐹𝑖,𝑗+1𝑎,𝑓,𝛾𝑓>0,𝑖,𝑗,𝑘,𝑙𝑤,𝑓,𝑒′,𝐹𝑖,𝑗+1𝑎,𝑓,𝛾<0.(3.37) The volumetric transfer across the matrix-fracture interface 𝑒 is given by𝑄𝑖,𝑗,𝑘,𝑙𝑤,𝑓,𝑒𝐹=−𝑖,𝑗+1𝑎,𝑚,𝑒,𝐾𝑓𝑖,𝑗,𝑘,𝑙𝑤,𝑒,𝐾+𝐹𝑖,𝑗+1𝑎,𝑚,𝑒,𝐾′𝑓𝑖,𝑗,𝑘,𝑙𝑤,𝑒,𝐾′𝜀,(3.38) where𝑓𝑖,𝑗,𝑘,𝑙𝑤,𝑒,𝐾=𝑓𝑖,𝑗,𝑘𝑤,𝑚,𝐾,𝐹𝑖,𝑗+1𝑎,𝑚,𝑒,𝐾𝑓>0,𝑖,𝑗,𝑘,𝑙𝑤,𝑓,𝑒,𝐹𝑖,𝑗+1𝑎,𝑚,𝑒,𝐾<0.(3.39) 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 matrix-fracture interfaces, and then obtain the matrix-vector form of (3.36)1Δ𝑡𝑖,𝑗,𝑘,𝑙𝐌𝑓𝐒𝑖,𝑗,𝑘,𝑙+1𝑤,𝑓−𝐒𝑖,𝑗,𝑘,𝑙𝑤,𝑓+𝐀𝑖,𝑗+1𝑠,𝑓,𝑓𝐟𝑖,𝑗,𝑘,𝑙𝑤,𝑓+𝐀𝑖,𝑗+1𝑠,𝑓,𝑚𝐟𝑖,𝑗,𝑘𝑤,𝑚=𝐐𝑖,𝑗,𝑘,𝑙𝑠,𝑓.(3.40) After 𝑁𝑠,𝑓 smaller time-steps, we update the saturation in the matrix domain by employing the following matrix-vector form of (3.32) as 1Δ𝑡𝑖,𝑗,𝑘𝐌𝑚𝐒𝑖,𝑗,𝑘+1𝑤,𝑚−𝐒𝑖,𝑗,𝑘𝑤,𝑚+𝐀𝑖,𝑗+1𝑠,𝑚,𝑚𝐟𝑖,𝑗,𝑘𝑤,𝑚+𝐀𝑖,𝑗+1𝑠,𝑚,𝑓̂𝐟𝑖,𝑗,𝑘𝑤,𝑓=𝐐𝑖,𝑗,𝑘𝑠,𝑚.(3.41)

In (3.40) and (3.41), 𝐀𝑖,𝑗+1𝑠,𝑚,𝑓 and 𝐀𝑖,𝑗+1𝑠,𝑓,𝑚 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 time-step 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𝑆𝑒=𝑆𝑤−𝑆𝑟𝑤1−𝑆𝑟𝑤−𝑆𝑟𝑛,(4.1) where 𝑆𝑟𝑤 and 𝑆𝑟𝑛 are the residual saturations of the wetting and nonwetting phases, respectively. The capillary pressure function [1] is given by𝑝𝑐𝑆𝑤=−𝐵𝑐𝑆log𝑒,(4.2) where 𝐵𝑐 is a positive parameter related to the absolute permeability. The relative permeabilities of two phases are given by𝑘𝑟𝑤=𝑆𝑑𝑒,𝑘𝑟𝑛=1−𝑆𝑒𝑑,(4.3) 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 right-hand 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 time-step 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.

Table 1 
Relevant data for Examples 1–4.
Table 2 
Comparison of Example 1.
Table 3 
Comparison of Example 2.
Table 4 
Comparison of Example 3.
Table 5 
Comparison of Example 4.
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 105 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.


Figure 1 
Distribution of fractures: Example 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 time-step 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.


Figure 2 
Water saturation profile at 0.35 PVI by the proposed method: Example 1.

Figure 3 
Water saturation profile at 0.35 PVI by the conventional method: Example 1.

In this example, the two-scale time splitting strategy is used only for the pressure equation, but with the single-scale time-stepping for the saturation equation. The two methods have the same subtiming between the pressure and the saturation, that is, 𝑁𝑠,𝑚=6. In order to achieve stable solutions, with two-scale strategy for the pressure equation, the proposed method can take the time-step size 0.9827 days, which is twice as many over the maximum time-step size 0.4121 days required by the conventional method. It is evident that the size of the pressure equation in the fine-scale time-steps is less than that in the coarse scale. Hence, this example demonstrates the efficiency of the proposed two-scale 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 106 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.


Figure 4 
Distribution of fractures: Example 2.

In this example, the proposed method takes the two-scale 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.


Figure 5 
Water saturation profile at 0.45 PVI by the proposed method: Example 2.

Figure 6 
Water saturation profile at 0.45 PVI by the conventional method: Example 2.

As shown in Table 3, in order to achieve stable solutions, the proposed method, with two-scale subtiming, can take a very large time-step size 1.6846 days, while the conventional method requires a very small time-step size, which is less than a tenth part of the proposed method. This example demonstrates that the two-scale 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.


Figure 7 
Distribution of fractures: Example 3.

In this example, we use the proposed method with two-scale 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.


Figure 8 
Water saturation profile at 0.5 PVI by the proposed method: Example 3.

Figure 9 
Water saturation profile at 0.5 PVI by the conventional method: Example 3.

From Table 4, we can see that the multiscale time splitting technique is capable of taking a very large time-step size 1.9010 days to achieve stable solutions. The conventional method, however, requires a very small time-step size, that is, 0.0895 days, which is the maximum time-step size for its stability with 𝑁𝑠,𝑚=5. Apparently, the multiscale time splitting strategy is superior to the conventional method for the multiple-fractured 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 matrix-fracture 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.


Figure 10 
Water saturation profile at 0.5 PVI by the proposed method: Example 4.

Figure 11 
Water saturation profile at 0.5 PVI by the conventional method: Example 4.

5. Conclusions

A multiscale time splitting strategy is introduced for simulating two-phase 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 matrix-fracture interface results from the continuity of capillary pressure. As a result, the conventional single-scale temporal schemes have shortage to handle the complexity of the multiphysics processes of two-phase 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 time-step size to achieve stable solutions, but a very small time-step 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.