© H. Gao et al., published by IFP Energies nouvelles, 2020

1 Introduction

Modeling of two-phase flows in porous media plays a crucial role in petroleum reservoir engineering; for instance, in the secondary oil recovery, water is injected to displace oil from reservoirs. In fractured media, there exist two distinct pore spaces: the fractures and the matrix. Compared to the matrix, the fractures usually have the higher porosity and permeability, but occupy very small volumes. Modeling of two-phase flows in fractured media is a challenging task in petroleum reservoir engineering and subsurface flow [1–3].

The classical two-phase flow models [2, 4–7] have been developed based on the single-phase flow model through introducing several new quantities, such as saturation, relative permeability, and capillary pressure. Relative permeability considerably effects on oil and water two-phase flows in porous media. Measurements and estimation formulations on relative permeabilities have been widely studied and successfully applied in petroleum reservoir engineering. Capillary effect [2, 3, 6, 8–12] is frequently viewed as the leading mechanism of oil recovery in fractured reservoirs, which results from the surface tension between immiscible fluids at the pore scales [13]. In fact, at the pore scales, the interfaces always exist between any binary fluids. Interactions between fluid molecules on the interfaces lead to surface tension, which can be calculated by the thermodynamically consistent diffuse interface models [13–15]. Recently, the energy concept has been taken into account in the conventional capillary formulations of two-phase flows in porous media (see [16] and the references therein). A phase-field model for fracture propagation in porous media filled with two-phase fluids has been proposed in [17]. However, the energy reconstructed from the classical model involves the ideal contribution of the wetting-phase saturation only, but the ideal fluid energy of both phases and attraction effect between two phases have not been taken into account, so the classical two-phase flow model is not completely thermodynamically consistent. Phase-field models for two-phase flows have been extensively studied in the literature, for instance [18–22]. Thermodynamically consistent phase-field-based diffuse interface models and simulation of two-phase flows in porous media have been reported in the literature, for instance [23–30], which introduce free energy potentials to characterize capillarity effect caused by the surface tension. The key ingredient of phase-field models of binary fluids is that a free energy potential is introduced to characterize phase behaviors, especially capillarity effect caused by the surface tension. In phase-field models, the general free energy is usually composed of two terms: the bulk free energy and gradient energy. The latter is introduced to account for the gradual phase transmission between binary bulk fluids, which results in clear, but very thin diffusive interfaces rather than sharp interfaces. For a realistic two-phase system, the interface thickness is usually less than ten nanometers [13–15, 31, 32]. However, there are no clear diffuse interfaces for binary immiscible fluids in porous media; actually, clear interfaces only exist at the microscopic pore scales. The phase saturations, as well as Darcy’s velocities, are the volumetric variables at the macroscopic Darcy’s scales. This means that the gradient energy has no physical basis for two-phase flows in porous media. Moreover, the saturation may be discontinuous in heterogeneous media and in fractured media [3, 6], and thus the gradient of saturation may be not well defined from a mathematical point of view. In addition, the effect of relative permeabilities is never taken into account in some phase-field models of two-phase flows in porous media. To fix the above defects, we will propose a new energy-based model of two-phase flows in porous media, which will be proved to be thermodynamically consistent and have great potentials in practical engineering applications.

For multi-phase flows in fractured media, several mathematical models have been developed in the literature [2, 3, 6, 33–40], such as the single-porosity model, the dual-porosity/dual-permeability model, and the discrete-fracture model. The single-porosity model treats the fractures as the matrix, and thus the excessive gridcells are needed due to the geometrical scale contrast between the matrix and the fractures. In the dual-porosity/dual-permeability model [34, 37, 39, 40], empirical functions to describe the matrix-fracture mass transfer incorporate ad hoc shape factors, and moreover, no theory is available for the calculation of shape factors in two-phase flows with capillary [3]. A two-phase flow model with local nonequilibrium in double porosity media has been derived in [41] using the homogenization approach, and this model includes the long-memory effects of the microscopic disequilibrium and the mass transfer between fractures and blocks. The discrete fracture model [35, 36] is based on the conception that the fracture aperture is very small compared to the matrix blocks, and as a consequence, we can simplify the fracture as the lower dimensional domain to reduce the contrast of geometric scales occurring in the single-porosity model. It is demonstrated in [3, 6] that the discrete fracture model is more preferable in computational efficiency. In this paper, we will present a discrete fracture energy-based model for two-phase flows in fractured media.

The IMplicit Pressure Explicit Saturation (IMPES) method [42–51] is an extensively used time-stepping approach for multi-phase flow simulation. Based on multi-physics processes of two-phase flows, the IMPES method treats the pressure implicitly and the saturation explicitly in the pressure equation, and thus leads to a linear pressure equation. The treatments result in the decoupling relationship between the pressure equation and the saturation equation, and as a result, the saturation can be calculated explicitly. In this paper, an IMPES scheme will be developed for the proposed model of two-phase flows in heterogeneous media and fractured media.

The rest of this paper is structured as follows. In Section 2, we present the model formulations for two-phase flows in heterogeneous and fractured media. In Section 3, the IMPES scheme is formulated for the proposed model. In Section 4, the experimental verification of the model is given. In Section 5, numerical examples are provided and discussed. Finally, the concluding remarks are provided in Section 6.

2 Mathematical model

In this section, we present a thermodynamically consistent energy-based model of incompressible, immiscible two-phase flows in heterogeneous porous media and fractured media.

2.1 Basic model equations

Here, the wetting phase is denoted by the subscript w, while the non-wetting phase is denoted by the subscript n. For the oil/water two-phase system, water and oil are the wetting phase and non-wetting phase, respectively. We let ϕ be the porosity of a porous medium, which may vary in space, but keeps constant with time. The saturation of phase α, denoted by S _α , where α = w, n, is the volume fraction of the void of a porous medium filled by phase α. The void volume of porous medium is jointly filled by binary fluids, so the saturations satisfy the constraint,

$$S_w+S_n=1.$$(1)

The mass conservation equation of each fluid is formulated as,

$$\varphi \frac{\mathrm{\partial }{S}_{\alpha }}{\mathrm{\partial }t}+\nabla \cdot {\mathbf{u}}_{\alpha }=q_{\alpha },\hspace{1em}\alpha =w,n,$$(2)where t is the time, u _α denotes Darcy’s velocity of phase α, and q _α stands for the injection/production rate of phase α.

The phase-field models introduce the free energy potentials to characterize phase behaviors, and capillarity effect caused by the surface tension is also modeled therein. In phase-field models, the general free energy usually contains the bulk free energy and gradient energy. The latter accounts for the gradient energy contributions on the diffuse interfaces. For a realistic two-phase system, the interface thickness is usually no more than a few nanometers [13–15, 32]. So the diffuse interfaces between two phases only exist at the microscopic pore scales, and the gradient energy has no physical meaning for two-phase flows in porous media. What’s more, the saturation may be discontinuous in heterogeneous and fractured media [3, 6], and thus the gradient of saturation may be not well defined from a mathematical point of view. Consequently, the gradient free energy in the classical phase-field models should not be taken into account for two-phase flows in porous media.

There have been several bulk free energy functions used for phase-field models in the literature, for instance, the double well potential and logarithmic Flory–Huggins energy potential. From a physical view of point, the logarithmic energy potential is more realistic than the double well potential [52–54]; in fact, it can be viewed as a simplified approximation of the Helmholtz free energy [55] determined by realistic equations of state, e.g. van der Waals equation of state and Peng–Robinson equation of state [56], which are widely employed in physics and oil–gas engineering respectively. Let S _rα be the residual saturation of phase α and let ς = 1−S _rw −S _rn . We propose the following logarithmic free energy function F as,

$$F\left(S_w,S_n\right)=\sum _{\alpha =w,n}\mathrm{}\varsigma {\gamma }_{\alpha }{S}_{\alpha }^e\left(\ln \left(S_{\alpha }^e\right)-1\right)+\varsigma {\gamma }_{\mathrm{wn}}{S}_w^e{S}_n^e,$$(3)where γ _α and γ _wn are the energy parameters that measure the capillary intensity, and $S_{\alpha }^e$ is the normalized saturation defined by,

$$S_{\alpha }^e=\frac{S_{\alpha }-S_{\mathrm{r\alpha }}}{\varsigma }.$$

The free energy in equation (3) is a combination of Flory–Huggins energy potential and Helmholtz free energy determined by Peng–Robinson equation of state. The first part in the energy function equation (3) accounts for the ideal fluid energy, while the attraction effect is taken into account in the second part of F [55]. The gravitational potential energy is expressed as,

$$G(S_w,S_n)=\mathrm{gz}\sum _{\alpha =w,n}{S}_{\alpha }{\rho }_{\alpha },$$(4)where g is the gravity acceleration, z is the depth, and ρ _α is the mass density of phase α. The total energy U is the sum of the free energy and gravitational potential energy,

$$U(S_w,S_n)=F(S_w,S_n)+G(S_w,S_n).$$(5)

We remark that γ _α and γ _wn may depend on the medium properties, such as porosity and permeability, and thus it may vary in different spaces of heterogeneous and fractured media. Indeed, if we take γ = γ _w  = γ _n  = γ _wn , γ can be viewed as the macroscopic mean value of microscopic capillary pressure caused by the surface tension between binary fluids. There have been a few formulations to describe the microscopic capillary pressure, for instance, the famous Young–Laplace equation:

$$\gamma =\frac{2\sigma }{\epsilon },$$(6)where σ is the surface tension and ε measures the pore radius at pore scales.

The chemical potentials are defined as the partial derivatives of F with respect to $S_w$ and $S_n$,

$${\mu }_w\left(S_w,S_n\right)=\frac{\mathrm{\partial }F\left(S_w,S_n\right)}{\mathrm{\partial }{S}_w}={\gamma }_w\ln \left(S_w^e\right)+{\gamma }_{\mathrm{wn}}{S}_n^e,$$(7)

$${\mu }_n\left(S_w,S_n\right)=\frac{\mathrm{\partial }F\left(S_w,S_n\right)}{\mathrm{\partial }{S}_n}={\gamma }_n\ln \left(S_n^e\right)+{\gamma }_{\mathrm{wn}}{S}_w^e.$$(8)

In heterogeneous media and fractured media, ${\mu }_{\alpha }$ may be continuous in space, whereas ${\gamma }_{\alpha }$ and ${\gamma }_{\mathrm{wn}}$ may take different values for different types of media, and as a consequence, the saturation may be discontinuous in space. The derivative of the gravitational potential energy $G$ with respect to $S_{\alpha }$ is expressed as,

$${\nu }_{\alpha }=\frac{\mathrm{\partial }G}{\mathrm{\partial }{S}_{\alpha }}={\rho }_{\alpha }\mathrm{gz}.$$(9)

The modified Darcy’s law for phase $\alpha$ is formulated as,

$${\mathbf{u}}_{\alpha }=-\frac{k_{\mathrm{r\alpha }}}{{\eta }_{\alpha }}\mathbf{K}\nabla (p+{\mu }_{\alpha }+{\nu }_{\alpha }),\hspace{1em}\alpha =w,n,$$(10)where $p$ is the pressure, $\mathbf{K}$ represents the absolute permeability tensor of the porous medium, $k_{\mathrm{r\alpha }}$ is the relative permeability of phase $\alpha$ and ${\eta }_{\alpha }$ is the viscosity of phase $\alpha$. We note that the physical meaning of the pressure $p$ is relevant to the free energy, which will be discussed in detail in Section 2.3. The absolute permeability $\mathbf{K}$ is the symmetric positive definite tensor. The relative permeability $k_{\mathrm{r\alpha }}$ of phase $\alpha$ is a nonnegative function of $S_{\alpha }$, for instance [3],

$$k_{\mathrm{rw}}=(S_w^e{)}^{l_w},\hspace{1em}{k}_{\mathrm{rn}}=(S_n^e{)}^{l_n},$$(11)where $l_w$ and $l_n$ are positive integers. It is distinguished from the traditional Darcy model of two-phase flows in porous media that the chemical potential is involved in equation (12) as one of primary driving forces in addition to pressure. To simplify the notations, we define the phase mobility ${\lambda }_{\alpha }=\frac{k_{\mathrm{r\alpha }}}{{\eta }_{\alpha }}$ and rewrite equation (10) as,

$${\mathbf{u}}_{\alpha }=-{\lambda }_{\alpha }\mathbf{K}\nabla (p+{\mu }_{\alpha }+{\nu }_{\alpha }),\hspace{1em}\alpha =w,n.$$(12)

In summary, the phase-field model of two-phase flows in porous media is a system of equations (2) and (12) associated with the constraint equation (1) and the formulations of the free energy and chemical potential. The proper boundary conditions are also required to complete the model, which rely on specific problems.

2.2 Thermodynamical consistency

The free energy of the proposed model has been constructed based on thermodynamical principles, and thus it is thermodynamically consistent. We now prove that the model obeys the second law of thermodynamics; that is, the total energy U of isothermal fluids in a closed domain would be decreasing with time until an equilibrium state is reached [55, 57]. More precisely, in a spatial domain Ω, we consider the simplified model,

$$\varphi \frac{\mathrm{\partial }{S}_{\alpha }}{\mathrm{\partial }t}+\nabla \cdot {\mathbf{u}}_{\alpha }=0,\hspace{1em}\alpha =w,n,$$(13)

$${\mathbf{u}}_{\alpha }=-{\lambda }_{\alpha }\mathbf{K}\nabla (p+{\mu }_{\alpha }+{\nu }_{\alpha }),\hspace{1em}\alpha =w,n,$$(14)associated with the boundary conditions,

$${\mathbf{u}}_{\alpha }\cdot \mathbf{n}=0,\hspace{1em}\mathrm{on}\hspace{1em}\mathrm{\partial }\mathrm{\Omega },\hspace{1em}\alpha =w,n,$$(15)where n is the outward unit normal vector to ∂Ω.

Using the constraint S _w  + S _n  = 1, we derive the divergence-free condition from equation (13) as,

$$\nabla \cdot {(\mathbf{u}}_w+{\mathbf{u}}_n)=0.$$(16)

We define the auxiliary phase velocity,

$${\mathbf{v}}_{\alpha }=-\mathbf{K}\nabla (p+{\mu }_{\alpha }+{\nu }_{\alpha }),\hspace{1em}\alpha =w,n.$$(17)

Then equation (14) is rewritten as,

$${\mathbf{u}}_{\alpha }={\lambda }_{\alpha }{\mathbf{v}}_{\alpha },\hspace{1em}\alpha =w,n.$$(18)

It is deduced from equations (17) and (18) using equation (16) that,

$$\begin{array}{cc}\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}{\lambda }_{\alpha }{\mathbf{v}}_{\alpha }\cdot {\mathbf{K}}^{-1}{\mathbf{v}}_{\alpha }\mathrm{d}x& =\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}{\mathbf{u}}_{\alpha }\cdot {\mathbf{K}}^{-1}{\mathbf{v}}_{\alpha }\mathrm{d}x\\ & =-\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}{\mathbf{u}}_{\alpha }\cdot \nabla (p+{\mu }_{\alpha }+{\nu }_{\alpha })\mathrm{d}x\\ \begin{array}{c}\\ \end{array}& \begin{array}{c}=\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}\left(p+{\mu }_{\alpha }+{\nu }_{\alpha }\right)\nabla \bullet {\mathbf{u}}_{\alpha }\mathrm{d}x\\ =\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}\left({\mu }_{\alpha }+{\nu }_{\alpha }\right)\nabla \bullet {\mathbf{u}}_{\alpha }\mathrm{d}x\end{array}\end{array}$$(19)

Since λ _α  ≥ 0 and K is the symmetric positive definite tensor, we deduce using equation (13) that,

$$\begin{array}{cc}{\int }_{\mathrm{\Omega }}\mathrm{}\varphi \frac{\mathrm{\partial }U}{\mathrm{\partial }t}\mathrm{d}x& ={\int }_{\mathrm{\Omega }}\mathrm{}\varphi \frac{\mathrm{\partial }(F+G)}{\mathrm{\partial }t}\mathrm{d}x\\ & =\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}\varphi \frac{\mathrm{\partial }(F+G)}{\mathrm{\partial }{S}_{\alpha }}\frac{\mathrm{\partial }{S}_{\alpha }}{\mathrm{\partial }t}\mathrm{d}x\\ \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}=\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}({\mu }_{\alpha }+{\nu }_{\alpha })\varphi \frac{\mathrm{\partial }{S}_{\alpha }}{\mathrm{\partial }t}\mathrm{d}x\\ =-\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}\left({\mu }_{\alpha }+{\nu }_{\alpha }\right)\nabla \cdot {\mathbf{u}}_{\alpha }\mathrm{d}x\\ =-\sum _{\alpha =w,n}\mathrm{}{\int }_{\mathrm{\Omega }}\mathrm{}{\lambda }_{\alpha }{\mathbf{v}}_{\alpha }\cdot {\mathbf{K}}^{-1}{\mathbf{v}}_{\alpha }\mathrm{dx}\le 0\end{array}\end{array}.$$(20)

Thus, the proposed model is thermodynamically consistent.

2.3 Special cases

We now present several special cases of the proposed model, one of which is the classical two-phase flow model. The α-phase pressure is denoted by p _α , where α = w, n. The capillary pressure p _c is defined as the difference between two phase pressures as,

$$p_c=p_n-p_w.$$(21)

As usual, the capillary pressure can be described as a function of saturations. For the proposed model, we take,

$$p_c={\mu }_n-{\mu }_w.$$(22)

In the first special case, we take γ _n  = γ _wn  = 0, and as a result, μ _n  = 0. From equation (22), the capillary pressure changes into $p_c=-{\mu }_w$. Darcy’s velocity of the wetting phase is rewritten as,

$$\begin{array}{cc}{\mathbf{u}}_w& =-{\lambda }_w\mathbf{K}\nabla (p+{\mu }_w+{\nu }_w)\\ & =-{\lambda }_w\mathbf{K}\nabla (p-p_c+{\nu }_w)\\ & =-{\lambda }_w\mathbf{K}\nabla (p+p_w-p_n+{\nu }_w),\end{array}$$(23)which leads to $p=p_n$ as a result of,

$${\mathbf{u}}_w=-{\lambda }_w\mathbf{K}\nabla (p_w+{\nu }_w).$$(24)

Thus, we obtain the Darcy velocity of the non-wetting phase as,

$${\mathbf{u}}_n=-{\lambda }_n\mathbf{K}\nabla (p+{\nu }_n)=-{\lambda }_n\mathbf{K}\nabla (p_n+{\nu }_n).$$(25)

From equation (7), the capillary pressure has the form $p_c=-{\gamma }_w\ln \left(S_w^e\right)$. So this case is the classical two-phase flow model [3, 6, 42, 46–48]. In this case, the model usually describes an imbibition displacement process; for instance, water is displacing oil.

In the second special case, we take ${\gamma }_w={\gamma }_{\mathrm{wn}}=0$, and then the capillary pressure becomes $p_c={\mu }_n$. Darcy’s velocity of the non-wetting phase is rewritten as,

$$\begin{array}{cc}{\mathbf{u}}_n& =-{\lambda }_n\mathbf{K}\nabla (p+{\mu }_n+{\nu }_n)\\ & =-{\lambda }_n\mathbf{K}\nabla (p+p_c+{\nu }_n)\\ & =-{\lambda }_n\mathbf{K}\nabla (p+p_n-p_w+{\nu }_n)\end{array},$$(26)which results in $p=p_w$. Consequently, we obtain the formulation of Darcy’s velocities as,

$${\mathbf{u}}_{\mathrm{\alpha }}=-{\lambda }_{\alpha }\mathbf{K}\nabla (p_{\alpha }+{\nu }_{\alpha }),\hspace{1em}\alpha =w,n.$$(27)

In this case, the model can describe a drainage displacement process, for instance, oil displacing water [2].

In the third special case, we take ${\gamma }_w={\gamma }_n=\gamma$ and ${\gamma }_{\mathrm{wn}}=0$, where $\gamma$ is a constant. Moreover, $S_{\mathrm{rw}}=S_{\mathrm{rn}}=0$. In this case, we have ${\mu }_{\alpha }=\gamma \ln \left(S_{\alpha }\right)$ and then we deduce,

$$S_w\nabla {\mu }_w+S_n\nabla {\mu }_n=\gamma S_w\nabla \ln \left(S_w\right)+\gamma S_n\nabla \ln \left(S_n\right)=\gamma \nabla S_w+\gamma \nabla S_n=0.$$(28)

Taking into account equations (21), (22) and (28), we deduce the formulations of Darcy’s velocities as,

$$\begin{array}{cc}{\mathbf{u}}_w& =-{\lambda }_w\mathbf{K}\nabla (p+{\mu }_w+{\nu }_w)\\ & =-{\lambda }_w\mathbf{K}(\nabla p+\nabla {\mu }_w-S_n\nabla {\mu }_n-S_w\nabla {\mu }_w+\nabla {\nu }_w)\\ \begin{array}{c}\\ \\ \begin{array}{c}\\ \end{array}\end{array}& \begin{array}{c}=-{\lambda }_w\mathbf{K}(\nabla p+S_n\nabla ({\mu }_w-{\mu }_n)+\nabla {\nu }_w)\\ =-{\lambda }_w\mathbf{K}(\nabla p-S_n\nabla p_c+\nabla {\nu }_w)\\ \begin{array}{c}=-{\lambda }_w\mathbf{K}(\nabla p-S_n\nabla (p_n-p_w)+\nabla {\nu }_w)\\ =-{\lambda }_w\mathbf{K}(\nabla p_w+\nabla p-S_n\nabla p_n-S_w\nabla p_w+\nabla {\nu }_w)\end{array}\end{array}\end{array},$$(29)

$$\begin{array}{cc}{\mathbf{u}}_n& =-{\lambda }_n\mathbf{K}\nabla (p+{\mu }_n+{\nu }_n)\\ & =-{\lambda }_n\mathbf{K}(\nabla p+\nabla {\mu }_n-S_n\nabla {\mu }_n-S_w\nabla {\mu }_w+\nabla {\nu }_n)\\ \begin{array}{c}\\ \\ \begin{array}{c}\\ \end{array}\end{array}& \begin{array}{c}=-{\lambda }_n\mathbf{K}(\nabla p+S_w\nabla ({\mu }_n-{\mu }_w)+\nabla {\nu }_n)\\ =-{\lambda }_n\mathbf{K}(\nabla p+S_w\nabla p_c+\nabla {\nu }_n)\\ \begin{array}{c}=-{\lambda }_n\mathbf{K}(\nabla p+S_w\nabla (p_n-p_w)+\nabla {\nu }_n)\\ =-{\lambda }_n\mathbf{K}(\nabla p_n+\nabla p-S_n\nabla p_n-S_w\nabla p_w+\nabla {\nu }_n)\end{array}\end{array}\end{array},$$(30)which lead to the weighted pressure as,

$$\nabla p=S_n\nabla p_n+S_w\nabla p_w,$$(31)as a result of,

$${\mathbf{u}}_{\alpha }=-{\lambda }_{\alpha }\mathbf{K}\nabla (p_{\alpha }+{\nu }_{\alpha }),\hspace{1em}\alpha =w,n.$$(32)

In this case, the model describes the imbibition and drainage mixed displacement process.

The above derivations of special cases show that specified pressure and capillary pressure can be derived from specified free energy. This means that for the energy-based model proposed in this paper, the major effort will be devoted to designing the physically reasonable free energy formulations rather than caring the specific pressures and capillary pressure. Generally speaking, the free energy can be constructed using relevant physical laws and relations rather than simple empirical formulations. The resultant model will be naturally consistent with thermodynamics. This is a pronounced feature of the proposed energy-based model.

2.4 Alternative formulation

There is not a pressure equation in the system of two-phase flow model, so we need to reconstruct an equation for pressure. Following the two-phase flow formulation for the classical two-phase flow model [3, 6], we provide an alternative formulation for the proposed model. The alternative formulation is a system composed of a pressure equation and a saturation equation, which can be solved more conveniently than the original system.

We define the general phase potential,

$${\mathrm{\Phi }}_{\alpha }=p+{\mu }_{\alpha }+{\nu }_{\alpha },\hspace{1em}\alpha =w,n,$$(33)and the general capillary potential,

$${\mathrm{\Phi }}_c={\mathrm{\Phi }}_n-{\mathrm{\Phi }}_w={\mu }_n-{\mu }_w+{\nu }_n-{\nu }_w.$$(34)

Then we rewrite the velocity of phase $\alpha$ as,

$${\mathbf{u}}_{\alpha }=-{\lambda }_{\alpha }\mathbf{K}\nabla {\mathrm{\Phi }}_{\alpha },\hspace{1em}\alpha =w,n.$$(35)

The total velocity $u_t$ is defined as the sum of velocities of two phases: $u_t=u_n+u_w$, which has the equivalent form,

$$\begin{array}{cc}{\mathbf{u}}_t& =-{\lambda }_w\mathbf{K}\nabla {\mathrm{\Phi }}_w-{\lambda }_n\mathbf{K}\nabla {\mathrm{\Phi }}_n\\ & =-{\lambda }_t\mathbf{K}\nabla {\mathrm{\Phi }}_w-{\lambda }_n\mathbf{K}\nabla ({\mathrm{\Phi }}_n-{\mathrm{\Phi }}_w)\\ & =-{\lambda }_t\mathbf{K}\nabla {\mathrm{\Phi }}_w-{\lambda }_n\mathbf{K}\nabla {\mathrm{\Phi }}_c\end{array},$$(36)where ${\lambda }_t={\lambda }_w+{\lambda }_n$ is the total mobility. Furthermore, we define,

$${\mathbf{u}}_a=-{\lambda }_t\mathbf{K}\nabla {\mathrm{\Phi }}_w,\hspace{1em}{\mathbf{u}}_c=-{\lambda }_n\mathbf{K}\nabla {\mathrm{\Phi }}_c.$$(36)

As a result, ${\mathbf{u}}_t$ is reformulated as,

$${\mathbf{u}}_t={\mathbf{u}}_a+{\mathbf{u}}_c.$$(38)

We define the fractional flow function ${\xi }_w={\lambda }_w/{\lambda }_t$. The wetting-phase velocity is rewritten as,

$${\mathbf{u}}_w={\xi }_w{\mathbf{u}}_a.$$(39)

The model is reformulated as the following governing equations,

$$\nabla \cdot {\mathbf{u}}_t=-\nabla \cdot {\lambda }_t\mathbf{K}\nabla {\mathrm{\Phi }}_w-\nabla \cdot {\lambda }_n\mathbf{K}\nabla {\mathrm{\Phi }}_c=q_t,$$(40)

$$\varphi \frac{\mathrm{\partial }{S}_w}{\mathrm{\partial }t}+\nabla \cdot \left({\xi }_w{\mathbf{u}}_a\right)=q_w,$$(41)where $q_t=q_w+q_n$.

The boundary $\mathrm{\partial }\mathrm{\Omega }$ of the spatial domain $\mathrm{\Omega }$ is decomposed into the Dirichlet part ${\mathrm{\Gamma }}_D$ and Neumann part ${\mathrm{\Gamma }}_N$, where $\mathrm{\partial }\mathrm{\Omega }={\mathrm{\Gamma }}_D\bigcup {\mathrm{\Gamma }}_N$. The boundary conditions associated to equations (40) and (41) are given by,

$${\mathrm{\Phi }}_w={\mathrm{\Phi }}_D\hspace{1em}\mathrm{on}\hspace{1em}{\mathrm{\Gamma }}_D,$$(42)

$$({\mathbf{u}}_a+{\mathbf{u}}_c)\cdot \mathbf{n}=q_N\hspace{1em}\mathrm{on}\hspace{1em}{\mathrm{\Gamma }}_N,$$(43)where $n$ is the outward unit normal vector to $\mathrm{\partial }\mathrm{\Omega }$. The wetting-phase saturation on the boundary $\mathrm{\partial }\mathrm{\Omega }$ is subject to,

$$S_w=S_w^N\hspace{1em}\mathrm{on}\hspace{1em}{\mathrm{\Gamma }}_N.$$(44)

The initial saturation is provided by,

$$S_w=S_w^0\hspace{1em}\mathrm{in}\hspace{1em}\mathrm{\Omega }.$$(45)

2.5 Discrete-fracture model

Here, the discrete-fracture model [3, 6] is applied to describe two-phase flows in fractured media based on the proposed model. The key idea of the discrete-fracture model [3, 6] is that for a n-dimensional matrix-fracture domain, the matrix is n-dimensional, and the fractures are approximated as the (n−1)-dimension.

We use the subscript m to denote the variables defined in the matrix; for example, ${\mathrm{\Phi }}_{w,m}$ denotes the form of ${\mathrm{\Phi }}_w$ in the matrix. The pressure equation (40) in the matrix changes into,

$$-\nabla \cdot {\lambda }_{t,m}{\mathbf{K}}_m\nabla {\mathrm{\Phi }}_{w,m}-\nabla \cdot {\lambda }_{n,m}{\mathbf{K}}_m\nabla {\mathrm{\Phi }}_{c,m}=q_{t,m}.$$(46)

In the fracture system, the potentials and saturations are constant along the fracture width, and the pressure equation in the fractures is expressed as,

$$-\nabla \cdot {\lambda }_{t,f}{\mathbf{K}}_f\nabla {\mathrm{\Phi }}_{w,f}-\nabla \cdot {\lambda }_{n,f}{\mathbf{K}}_f\nabla {\mathrm{\Phi }}_{c,f}=q_{t,f}+Q_{t,f},$$(47)where the subscript f denotes the fracture and $Q_{t,f}$ stands for the total mass transfer across the matrix-fracture boundaries. The velocities ${\mathbf{u}}_a$ in the matrix and fractures are expressed as,

$${\mathbf{u}}_{a,m}=-{\lambda }_{t,m}{\mathbf{K}}_m\nabla {\mathrm{\Phi }}_{w,m},\hspace{1em}{\mathbf{u}}_{a,f}=-{\lambda }_{t,f}{\mathbf{K}}_f\nabla {\mathrm{\Phi }}_{w,f}.$$(48)

Similarly, we can formulate the saturation equation in the matrix as,

$${\varphi }_m\frac{\mathrm{\partial }{S}_{w,m}}{\mathrm{\partial }t}+\nabla \cdot \left({\xi }_{w,m}{\mathbf{u}}_{a,m}\right)=q_{w,m},$$(49)and the saturation equation in the fractures as,

$${\varphi }_f\frac{\mathrm{\partial }{S}_{w,f}}{\mathrm{\partial }t}+\nabla \cdot \left({\xi }_{w,f}{\mathbf{u}}_{a,f}\right)=q_{w,f}+Q_{w,f},$$(50)where $Q_{w,f}$ represents the mass transfer of wetting phase across the matrix-fracture boundaries.

3 Numerical method

The IMplicit Pressure Explicit Saturation (IMPES) method [2, 42–44, 50] is a splitting approach based on physics, in which only a pressure equation is solved implicitly, and the saturation is updated explicitly as long as the pressure is obtained. Compared to the fully implicit scheme, the IMPES method requires smaller computational cost and memory at each time step. This advantage is more pronounced for large-scale numerical simulations.

The total time interval $[0,T]$ is divided into $M$ time steps as $0=t_0<t_1<\dots <t_M=T$. We denote the time step length by $\Delta t_i=t_{i+1}-t_i$. The value of a function $v$ at the time $t_i$ is denoted by $v^i$.

3.1 General IMPES method

Following the IMPES method for the classical two-phase flow model [3, 6], we present the IMPES scheme for the proposed model. The relative permeabilities $k_{\mathrm{r\alpha }}$ are calculated from the saturations at the previous time step, and thus the variables ${\lambda }_w$, ${\lambda }_n$, and ${\lambda }_t$ in the pressure equation are treated explicitly. The chemical potentials are computed explicitly as well. The pressure equation is semi-implicitly discretized as,

$$-\nabla \cdot {\lambda }_t^i\mathbf{K}\nabla {\mathrm{\Phi }}_w^{i+1}-\nabla \cdot {\lambda }_n^i\mathbf{K}\nabla {\mathrm{\Phi }}_c^i=q_t^{i+1},$$(51)which is a linear, symmetric positive definite equation of pressure and can be easy to be solved. Once the pressure is calculated, the velocity ${\mathbf{u}}_a^{i+1}$ can be evaluated explicitly as,

$${\mathbf{u}}_a^{i+1}=-{\lambda }_t^i\mathbf{K}\nabla {\mathrm{\Phi }}_w^{i+1}.$$(52)

The saturation equation of the wetting phase is discretized by the forward Euler time scheme as,

$$\varphi \frac{S_w^{i+1}-S_w^i}{\Delta t_i}+\nabla \cdot \left({\xi }_w^i{\mathbf{u}}_a^{i+1}\right)=q_w^{i+1}.$$(53)

As a result, the saturations are updated explicitly.

3.2 IMPES method for the discrete-fracture model

The semi-implicit time schemes for the pressure equations (46) and (47) in the matrix domain and fractures are formulated as,

$$-\nabla \cdot {\lambda }_{t,m}^i{\mathbf{K}}_m\nabla {\mathrm{\Phi }}_{w,m}^{i+1}-\nabla \cdot {\lambda }_{n,m}^i{\mathbf{K}}_m\nabla {\mathrm{\Phi }}_{c,m}^i=q_{t,m}^{i+1},$$(54)

$$-\nabla \cdot {\lambda }_{t,f}^i{\mathbf{K}}_f\nabla p_f^{i+1}-\nabla \cdot {\lambda }_{n,f}^i{\mathbf{K}}_f\nabla {\mathrm{\Phi }}_{c,f}^i=q_{t,f}^{i+1}+Q_{t,f}^{i+1}.$$(55)

Once the pressure equations are solved, the velocities are calculated as,

$${\mathbf{u}}_{a,m}^{i+1}=-{\lambda }_{t,m}^i{\mathbf{K}}_m\nabla {\mathrm{\Phi }}_{w,m}^{i+1},\hspace{1em}{\mathbf{u}}_{a,f}^{i+1}=-{\lambda }_{t,f}^i{\mathbf{K}}_f\nabla {\mathrm{\Phi }}_{w,f}^{i+1}.$$(56)

The saturations in the matrix domain and fractures are updated as,

$${\varphi }_m\frac{S_{w,m}^{i+1}-S_{w,m}^i}{\Delta t_i}+\nabla \cdot \left({\xi }_{w,m}^i{\mathbf{u}}_{a,m}^{i+1}\right)=q_{w,m}^{i+1},$$(57)

$${\varphi }_f\frac{S_{w,f}^{i+1}-S_{w,f}^i}{\Delta t_i}+\nabla \cdot \left({\xi }_{w,f}^i{\mathbf{u}}_{a,f}^{i+1}\right)=q_{w,f}^{i+1}+Q_{w,f}^{i+1}.$$(58)

4 Model verification

In this section, we verify the proposed model using the experimental data. We recall that a physically reasonable free energy has already introduced in this paper, and then the corresponding chemical potentials are included in the modified Darcy law. As shown in equation (22), the capillary pressure can be obtained from the difference between chemical potentials of two phases. To validate the model, we employ the capillary pressure formulation equation (22) to fit the experimental data of water–oil systems provided in [58]. Using the constraint $S_w+S_n=1$, we can express the capillary pressure function given in equation (22) as a function of $S_w$, i.e., $p_c\left(S_w\right)$. The energy parameters can be calculated by the least squares method: solve the minimizer of,

$$J\left({\gamma }_w,{\gamma }_n,{\gamma }_{\mathrm{nw}}\right)=\sum _{i=1}^m\mathrm{}{\left(p_c\left(S_{w,i}\right)-{\bar{p}}_{c,i}\right)}^2,$$(59)where $m$ is the sample size and ${\bar{p}}_{c,i}$ is the experimental value of capillary pressure measured at $S_{w,i}$.

Figure 1 shows the capability and performance of the capillary pressure function equation (22) in approximating the experimental data of three samples in [58]. Here, the wetting-phase residual saturations are taken as $S_{\mathrm{rw}}=0.098, 0.039, 0.03$ for three samples respectively, while $S_{\mathrm{rn}}=0$ is taken for all samples. It can be observed that the capillary pressure functions are essentially in agreement with the experimental data. The capillary pressure curves clearly demonstrate that the capillary pressure functions are non-convex functions with respect to the saturation $S_w$, and thus the conventional convex functions [2, 4–7, 16] for the capillary pressure would not agree with the experimental data.

5 Numerical results

The Cell-Centered Finite Difference (CCFD) method is used as the spatial discretization method in numerical simulation. The detailed formulations of CCFD for the discrete-fracture model can be found in [59].

In numerical examples, we consider the horizontal layers in fractured media with the spatial domain dimensions 10 m × 10 m × 1 m, but containing different fracture networks. The fracture aperture is 0.01 m. The porosities of the matrix and fractures are taken as ${\varphi }_m=0.2$ and ${\varphi }_f=1$. The absolute permeabilities in the matrix blocks and the fractures are 1 md and 10⁵ md, respectively. The quadratic relative permeabilities are applied in both fractures and matrix; that is, $l_w=l_n=2$ in equation (11). The viscosities of water and oil are ${\eta }_w=1$ cP and ${\eta }_n=0.5$ cP, respectively. The residual saturations are taken as $S_{\mathrm{rw}}=0.03$ and $S_{\mathrm{rn}}=0$. The energy parameters are taken as ${\gamma }_w^m=0.1166$ bar, ${\gamma }_n^m=0.0108$ bar, ${\gamma }_{\mathrm{nw}}^m=0.0742$ bar in the matrix and ${\gamma }_w^f=1.166\mathrm{e}-3$ bar, ${\gamma }_n^f=1.08\mathrm{e}-4$ bar, ${\gamma }_{\mathrm{nw}}^f=7.42\mathrm{e}-4$ bar in the fractures. The effect of gravity is neglected in the horizontal media.

5.1 Example 1

In this example, we verify thermodynamical consistency of the proposed model; that is, the total free energy of isothermal fluids in a closed domain is decreasing with time. The mass fluxes towards outsides of all boundaries vanish. There is no injection or extraction to the interior of the domain. The distribution of three parallel fractures in the fractured medium is illustrated in Figure 2.

The diffusion between oil and water takes place due to the capillarity effect. The dynamical process is simulated until 1 Pore Volume Injection (PVI). The water distributions at different times are shown in Figure 3. The variation of total energy with time is depicted in Figure 4. It is apparent that total energy is decreasing with time, and as a consequence, thermodynamical consistency of the proposed model is verified.

5.2 Example 2

The fractured medium is the same to that in Example 1 as illustrated in Figure 2. The void of the media is initially saturated with oil and residual water ($S_{\mathrm{rw}}=0.03$). Water is injected at the left end of the medium, while oil is produced at the right-hand side. The injection rate is 0.2 PV/year. The mass fluxes towards outsides of the other boundaries vanish. There is no other injection or extraction to the interior of the domain.

The displacement process of oil by water is simulated until 0.39 PVI. The water distributions at different times are shown in Figure 5. For comparison, the case without free energy, i.e., with zero energy parameters, is simulated as well and the corresponding results are depicted in Figure 6.

Figure 6 shows that water flows quickly through the fractures and the breakthrough takes place early. In contrast, Figure 5 demonstrates that the capillarity effect caused by surface tension significantly delays breakthrough, and thus it leads to more efficiency of the displacement of oil by water. Figure 7 illustrates that the free energy in the regions occupied by oil and fractures is clearly higher than the other regions, which leads to the energy transmission and dissipation, and as a result, water tends to flow from fractures to matrix regions. This reveals that capillarity effect is a result of the free energy transmission and dissipation.

5.3 Example 3

In this example, we consider a fractured medium with a fracture network that is illustrated in Figure 8. The displacement process of oil by water is simulated until 0.33 PVI. The initial and boundary conditions are the same to Example 2. The water distributions at different times are depicted in Figure 9. For comparison, the results without free energy are also shown in Figure 10. The energy distributions at different times are illustrated in Figure 11.

From Figure 9, we can see that in the displacement process, water is driven to cross from fractures to matrix regions due to the capillarity effect caused by surface tension. The contrast phenomena is observed in Figure 10 that water flows quickly through the fractures and it is accumulating at the ends of fractures. Consequently, capillarity can improve the oil recovery largely. Figure 11 illustrates that the fluids in the oil regions and fractures have the higher free energy, which leads to more water flowing from fractures to the matrix.

6 Conclusion

A thermodynamically consistent energy-based model has been developed for two-phase flows in heterogeneous and fractured media. The model inherits the framework of the traditional model of two-phase flows in porous media, such as relative permeability, volumetric phase velocity, and capillarity effect. A logarithmic energy potential, which is more realistic than the commonly used double well potential, is proposed to characterize the capillarity effect. We develop the energy-based discrete fracture model to describe two-phase flows in fractured media. The popularly used implicit pressure explicit saturation method is applied to simulate the model equations. The model is verified by the experimental data. Numerical simulation results are provided to demonstrate that the proposed model can provide physically reasonable results.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No. 41907167), the Natural Science Foundation of Zhejiang Province (No. LY18E090004), the Scientific and Technical Research Project of Hubei Provincial Department of Education (No. D20192703) and the Technology Creative Project of Excellent Middle & Young Team of Hubei Province (No. T201920).

References

1. Cai J., Wei W., Hu X., Liu R., Wang J. (2017) Fractal characterization of dynamic fracture network extension in porous media, Fractals 25, 1750023.
2. Chen Z., Huan G., Ma Y. (2006) Computational methods for multiphase flows in porous media, Society for Industrial and Applied Mathematics, Philadelphia, PA.
3. Hoteit H., Firoozabadi A. (2008) An efficient numerical model for incompressible two-phase flow in fractured media, Adv. Water Resour. 31, 891–905.
4. Chen J., Sun S., Wang X.-P. (2014) A numerical method for a model of two-phase flow in a coupled free flow and porous media system, J. Comput. Phys. 268, 1–16.
5. Galusinski C., Saad M. (2008) Two compressible immiscible fluids in porous media, J. Differ. Equ. 244, 1741–1783.
6. Hoteit H., Firoozabadi A. (2008) Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures, Adv. Water Resour. 31, 56–73.
7. Hou J., Chen J., Sun S., Chen Z. (2016) Adaptive mixed-hybrid and penalty discontinuous Galerkin method for two-phase flow in heterogeneous media, J. Comput. Appl. Math. 307, 262–283.
8. Cai J., Yu B. (2011) A discussion of the effect of tortuosity on the capillary imbibition in porous media, Transp. Porous Med. 89, 251–263.
9. Cai J., Yu B., Mei M., Luo L. (2010) Capillary rise in a single tortuous capillary, Chin. Phys. Lett. 27, 054701.
10. Shen A., Liu Y., Farouq Ali S.M. (2020) A model of spontaneous flow driven by capillary pressure in nanoporous media, Capillarity 3, 1–7.
11. Sun S. (2019) Darcy-scale phase equilibrium modeling with gravity and capillarity, J. Comput. Phys. 399, 108908.
12. Wei W., Cai J., Xiao J., Meng Q., Xiao B., Han Q. (2018) Kozeny-Carman constant of porous media: Insights from fractal-capillary imbibition theory, Fuel 234, 1373–1379.
13. Kou J., Sun S., Wang X. (2015) Efficient numerical methods for simulating surface tension of multi-component mixtures with the gradient theory of fluid interfaces, Comput. Methods Appl. Mech. Eng. 292, 92–106.
14. Miqueu C., Mendiboure B., Graciaa A., Lachaise J. (2003) Modelling of the surface tension of pure components with the gradient theory of fluid interfaces: a simple and accurate expression for the influence parameters, Fluid Phase Equilib. 207, 225–246.
15. Miqueu C., Mendiboure B., Graciaa A., Lachaise J. (2005) Modeling of the surface tension of multicomponent mixtures with the gradient theory of fluid interfaces, Ind. Eng. Chem. Res. 44, 3321–3329.
16. Cancès C. (2018) Energy stable numerical methods for porous media flow type problems, Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles 73, 78.
17. Lee S., Mikelić A., Wheeler M.F., Wick T. (2018) Phase-field modeling of two phase fluid filled fractures in a poroelastic medium, Multiscale Model. Simul. 16, 1542–1580.
18. Abels H., Garcke H., Grun G. (2012) Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities, Math. Mod. Meth. Appl. Sci. 22, 1150013.
19. Guo Z., Lin P. (2015) A thermodynamically consistent phase-field model for two-phase flows with thermocapillary effects, J. Fluid Mech. 766, 226–271.
20. Shen J., Yang X. (2014) Decoupled energy stable schemes for phase-field models of two-phase complex fluids, SIAM, J. Sci. Comput. 36, B122–B145.
21. Zhu G., Kou J., Yao B., Wu Y.-S., Yao J., Sun S. (2019) Thermodynamically consistent modelling of two-phase flows with moving contact line and soluble surfactants, J. Fluid Mech. 879, 327–359.
22. Zhu G., Li A. (2020) Interfacial dynamics with soluble surfactants: A phase-field two-phase flow model with variable densities, Adv. Geo-Energy Res. 4, 86–98.
23. Chen C.-Y., Yan P.-Y. (2015) A diffuse interface approach to injection-driven flow of different miscibility in heterogeneous porous media, Phys. Fluids 27, 083101.
24. Cogswell D.A., Szulczewski M.L. (2017) Simulation of incompressible two-phase flow in porous media with large timesteps, J. Comput. Phys. 345, 856–865.
25. Cueto-Felgueroso L., Juanes R. (2014) A phase-field model of two-phase Hele-Shaw flow, J. Fluid Mech. 758, 522–552.
26. Dede L., Garcke H., Lam K.F. (2018) A Hele-Shaw-Cahn-Hilliard model for incompressible two-phase flows with different densities, J. Math. Fluid Mech. 20, 531–567.
27. Feng X.B., Wise S. (2012) Analysis of a Darcy-Cahn-Hilliard diffuse interface model for the Hele-Shaw flow and its fully discrete finite element approximation, SIAM J. Numer. Anal. 50, 1320–1343.
28. Ngamsaad W., Yojina J., Triampo W. (2010) Theoretical studies of phase-separation kinetics in a Brinkman porous medium, J. Phys. A: Math. Theor. 43, 202001.
29. Sabooniha E., Rokhforouz M.R., Ayatollahi S. (2019) Pore-scale investigation of selective plugging mechanism in immiscible two-phase flow using phase-field method, Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles 74, 78.
30. Wise S.M. (2010) Unconditionally stable finite difference, nonlinear multigrid simulation of the Cahn-Hilliard-Hele-Shaw system of equations, J. Sci. Comput. 44, 38–68.
31. Chen J., Chung E.T., He Z., Sun S. (2020) Generalized multiscale approximation of mixed finite elements with velocity elimination for subsurface flow, J. Comput. Phys. 404, 109133.
32. Kou J., Sun S., Wang X. (2020) A novel energy factorization approach for the diffuse-interface model with Peng-Robinson equation of state, SIAM J. Sci. Comput. 42, B30–B56.
33. Ghorayeb K., Firoozabadi A. (2000) Numerical study of natural convection and diffusion in fractured porous media, SPE J. 5, 12–20.
34. Kazemi H., Merrill L.S. (1979) Numerical simulation of water imbibition in fractured cores, Old SPE J. 19, 175–182.
35. Lee S., Jensen C., Lough M. (2000) Efficient finite-difference model for flow in a reservoir with multiple length-scale fractures, SPE J. 5, 268–275.
36. Noorishad J., Mehran M. (1982) An upstream finite element method for solution of transient transport equation in fractured porous media, Water Resour. Res. 18, 588–596.
37. Pruess K. (1985) A practical method for modeling fluid and heat flow in fractured porous media, Old SPE J. 25, 14–26.
38. Sarkar S., Toksoz M.N., Burns D.R. (2002) Fluid flow simulation in fractured reservoirs, in: Report, Annual Consortium Meeting.
39. Thomas L., Dixon T., Pierson R. (1983) Fractured reservoir simulation, Old SPE J. 23, 42–54.
40. Warren J.E., Root P.J. (1963) The behavior of naturally fractured reservoirs, Old SPE J. 3, 245–255.
41. Amaziane B., Panfilov M., Pankratov L. (2016) Homogenized model of two-phase flow with local nonequilibrium in double porosity media, Adv. Math. Phys. 2016, 3058710.
42. Chen H., Kou J., Sun S., Zhang T. (2019) Fully mass-conservative IMPES schemes for incompressible two-phase flow in porous media, Comput. Methods Appl. Mech. Eng. 350, 641–663.
43. Chen Z., Huan G., Li B. (2004) An improved IMPES method for two-phase flow in porous media, Transp. Porous Med. 54, 361–376.
44. Coats K.H. (2003) IMPES stability: Selection of stable timesteps, SPE J. 8, 181–187.
45. Fagin R.G., Stewart C.H. Jr. (1966) A new approach to the two-dimensional multiphase reservoir simulator, Old SPE J. 6, 175–182.
46. Kou J., Sun S. (2010) A new treatment of capillarity to improve the stability of IMPES two-phase flow formulation, Comput. Fluids 39, 1293–1931.
47. Kou J., Sun S. (2013) Convergence of discontinuous Galerkin methods for incompressible two-phase flow in heterogeneous media, SIAM J. Numer. Anal. 51, 3280–3306.
48. Kou J., Sun S. (2014) Upwind discontinuous Galerkin methods with conservation of mass of both phases for incompressible two-phase flow in porous media, Numer. Methods Partial Differ. Equ. 30, 1674–1699.
49. Sheldon J.W., Zondek B., Cardwell W.T. (1959) One-dimensional, incompressible, noncapillary, two-phase fluid flow in a porous medium, Trans. SPE AIME 216, 290–296.
50. Yang H., Sun S., Yang C. (2017) Nonlinearly preconditioned semismooth Newton methods for variational inequality solution of two-phase flow in porous media, J. Comput. Phys. 332, 1–20.
51. Young L., Stephenson R. (1983) A generalized compositional approach for reservoir simulation, Old SPE J. 23, 727–742.
52. Cahn J.W., Allen S.M. (1977) A microscopic theory for domain wall motion and its experimental varification in fe-al alloy domain growth kinetics, J. Phys. Colloq. C7, C7–51.
53. Chen W., Wang C., Wang X., Wise S.M. (2019) Positivity-preserving, energy stable numerical schemes for the Cahn-Hilliard equation with logarithmic potential, J. Comput. Phys.: X 3, 100031.
54. Wang X., Kou J., Cai J. (2020) Stabilized energy factorization approach for Allen-Cahn equation with logarithmic Flory-Huggins potential, J. Sci. Comput. 82, 25.
55. Kou J., Sun S. (2018) Thermodynamically consistent modeling and simulation of multi-component two-phase flow with partial miscibility, Comput. Methods Appl. Mech. Eng. 331, 623–649.
56. Peng D., Robinson D.B. (1976) A new two-constant equation of state, Ind. Eng. Chem. Fund. 15, 59–64.
57. Kou J., Sun S. (2017) Efficient energy-stable dynamic modeling of compositional grading, Int. J. Numer. Anal. Mod. 14, 218–242.
58. Zhang C., Zhang C.-M., Zhang Z., Huang C., Mao Z. (2014) Effect of clay bound water on the mercury injection capillary pressure curves and the correction, Sci. Technol. Rev. 32, 44–49. (in Chinese).
59. Kou J., Sun S., Yu B. (2011) Multiscale time-splitting strategy for multiscale multiphysics processes of two-phase flow in fractured media, J. Appl. Math. 2011, 861905.

All Figures