Regular Article
Thermodynamically consistent modeling of twophase incompressible flows in heterogeneous and fractured media
^{1}
School of Civil Engineering, Shaoxing University, Shaoxing, 312000 Zhejiang, China
^{2}
Key Laboratory of Rock Mechanics and Geohazards of Zhejiang Province, Shaoxing, 312000 Zhejiang, China
^{3}
School of Mathematics and Statistics, Hubei Engineering University, Xiaogan, 432000 Hubei, China
^{4}
Computational Transport Phenomena Laboratory, Division of Physical Science and Engineering, King Abdullah University of Science and Technology, 239556900 Thuwal, Kingdom of Saudi Arabia
^{*} Corresponding authors: jishengkou@163.com; shuyu.sun@kaust.edu.sa
Received:
29
January
2020
Accepted:
2
April
2020
Numerical modeling of twophase flows in heterogeneous and fractured media is of great interest in petroleum reservoir engineering. The classical model for twophase flows in porous media is not completely thermodynamically consistent since the energy reconstructed from the capillary pressure does not involve the ideal fluid energy of both phases and attraction effect between two phases. On the other hand, the saturation may be discontinuous in heterogeneous and fractured media, and thus the saturation gradient may be not well defined. Consequently, the classical phasefield models can not be applied due to the use of diffuse interfaces. In this paper, we propose a new thermodynamically consistent energybased model for twophase flows in heterogeneous and fractured media, which is free of the gradient energy. Meanwhile, the model inherits the key features of the traditional models of twophase flows in porous media, including relative permeability, volumetric phase velocity and capillarity effect. To characterize the capillarity effect, a logarithmic energy potential is proposed as the free energy function, which is more realistic than the commonly used double well potential. The model combines with the discrete fracture model to describe twophase flows in fractured media. The popularly used implicit pressure explicit saturation method is used to simulate the model. Finally, the experimental verification of the model and numerical simulation results are provided.
© H. Gao et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Modeling of twophase 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 twophase flows in fractured media is a challenging task in petroleum reservoir engineering and subsurface flow [1–3].
The classical twophase flow models [2, 4–7] have been developed based on the singlephase flow model through introducing several new quantities, such as saturation, relative permeability, and capillary pressure. Relative permeability considerably effects on oil and water twophase 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 twophase flows in porous media (see [16] and the references therein). A phasefield model for fracture propagation in porous media filled with twophase fluids has been proposed in [17]. However, the energy reconstructed from the classical model involves the ideal contribution of the wettingphase 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 twophase flow model is not completely thermodynamically consistent. Phasefield models for twophase flows have been extensively studied in the literature, for instance [18–22]. Thermodynamically consistent phasefieldbased diffuse interface models and simulation of twophase 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 phasefield 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 phasefield 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 twophase 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 twophase 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 phasefield models of twophase flows in porous media. To fix the above defects, we will propose a new energybased model of twophase flows in porous media, which will be proved to be thermodynamically consistent and have great potentials in practical engineering applications.
For multiphase flows in fractured media, several mathematical models have been developed in the literature [2, 3, 6, 33–40], such as the singleporosity model, the dualporosity/dualpermeability model, and the discretefracture model. The singleporosity 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 dualporosity/dualpermeability model [34, 37, 39, 40], empirical functions to describe the matrixfracture mass transfer incorporate ad hoc shape factors, and moreover, no theory is available for the calculation of shape factors in twophase flows with capillary [3]. A twophase flow model with local nonequilibrium in double porosity media has been derived in [41] using the homogenization approach, and this model includes the longmemory 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 singleporosity 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 energybased model for twophase flows in fractured media.
The IMplicit Pressure Explicit Saturation (IMPES) method [42–51] is an extensively used timestepping approach for multiphase flow simulation. Based on multiphysics processes of twophase 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 twophase 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 twophase 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 energybased model of incompressible, immiscible twophase 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 nonwetting phase is denoted by the subscript n. For the oil/water twophase system, water and oil are the wetting phase and nonwetting 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,
The mass conservation equation of each fluid is formulated as,
(2)where t is the time, u _{ α } denotes Darcy’s velocity of phase α, and q _{ α } stands for the injection/production rate of phase α.
The phasefield models introduce the free energy potentials to characterize phase behaviors, and capillarity effect caused by the surface tension is also modeled therein. In phasefield 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 twophase 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 twophase 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 phasefield models should not be taken into account for twophase flows in porous media.
There have been several bulk free energy functions used for phasefield 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,
(3)where γ _{ α } and γ _{ wn } are the energy parameters that measure the capillary intensity, and is the normalized saturation defined by,
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,
(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,
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:
(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 and ,
In heterogeneous media and fractured media, may be continuous in space, whereas and 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 with respect to is expressed as,
The modified Darcy’s law for phase is formulated as,
(10)where is the pressure, represents the absolute permeability tensor of the porous medium, is the relative permeability of phase and is the viscosity of phase . We note that the physical meaning of the pressure is relevant to the free energy, which will be discussed in detail in Section 2.3. The absolute permeability is the symmetric positive definite tensor. The relative permeability of phase is a nonnegative function of , for instance [3],
(11)where and are positive integers. It is distinguished from the traditional Darcy model of twophase 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 and rewrite equation (10) as,
In summary, the phasefield model of twophase 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,
(14)associated with the boundary conditions,
(15)where n is the outward unit normal vector to ∂Ω.
Using the constraint S _{ w } + S _{ n } = 1, we derive the divergencefree condition from equation (13) as,
We define the auxiliary phase velocity,
Then equation (14) is rewritten as,
It is deduced from equations (17) and (18) using equation (16) that,
Since λ _{ α } ≥ 0 and K is the symmetric positive definite tensor, we deduce using equation (13) that,
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 twophase 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,
As usual, the capillary pressure can be described as a function of saturations. For the proposed model, we take,
In the first special case, we take γ _{ n } = γ _{ wn } = 0, and as a result, μ _{ n } = 0. From equation (22), the capillary pressure changes into . Darcy’s velocity of the wetting phase is rewritten as,
(23)which leads to as a result of,
Thus, we obtain the Darcy velocity of the nonwetting phase as,
From equation (7), the capillary pressure has the form . So this case is the classical twophase 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 , and then the capillary pressure becomes . Darcy’s velocity of the nonwetting phase is rewritten as,
(26)which results in . Consequently, we obtain the formulation of Darcy’s velocities as,
In this case, the model can describe a drainage displacement process, for instance, oil displacing water [2].
In the third special case, we take and , where is a constant. Moreover, . In this case, we have and then we deduce,
Taking into account equations (21), (22) and (28), we deduce the formulations of Darcy’s velocities as,
(30)which lead to the weighted pressure as,
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 energybased 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 energybased model.
2.4 Alternative formulation
There is not a pressure equation in the system of twophase flow model, so we need to reconstruct an equation for pressure. Following the twophase flow formulation for the classical twophase 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,
(33)and the general capillary potential,
Then we rewrite the velocity of phase as,
The total velocity is defined as the sum of velocities of two phases: , which has the equivalent form,
(36)where is the total mobility. Furthermore, we define,
As a result, is reformulated as,
We define the fractional flow function . The wettingphase velocity is rewritten as,
The model is reformulated as the following governing equations,
The boundary of the spatial domain is decomposed into the Dirichlet part and Neumann part , where . The boundary conditions associated to equations (40) and (41) are given by,
(43)where is the outward unit normal vector to . The wettingphase saturation on the boundary is subject to,
The initial saturation is provided by,
2.5 Discretefracture model
Here, the discretefracture model [3, 6] is applied to describe twophase flows in fractured media based on the proposed model. The key idea of the discretefracture model [3, 6] is that for a ndimensional matrixfracture domain, the matrix is ndimensional, 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, denotes the form of in the matrix. The pressure equation (40) in the matrix changes into,
In the fracture system, the potentials and saturations are constant along the fracture width, and the pressure equation in the fractures is expressed as,
(47)where the subscript f denotes the fracture and stands for the total mass transfer across the matrixfracture boundaries. The velocities in the matrix and fractures are expressed as,
Similarly, we can formulate the saturation equation in the matrix as,
(49)and the saturation equation in the fractures as,
(50)where represents the mass transfer of wetting phase across the matrixfracture 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 largescale numerical simulations.
The total time interval is divided into time steps as . We denote the time step length by . The value of a function at the time is denoted by .
3.1 General IMPES method
Following the IMPES method for the classical twophase flow model [3, 6], we present the IMPES scheme for the proposed model. The relative permeabilities are calculated from the saturations at the previous time step, and thus the variables , , and in the pressure equation are treated explicitly. The chemical potentials are computed explicitly as well. The pressure equation is semiimplicitly discretized as,
(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 can be evaluated explicitly as,
The saturation equation of the wetting phase is discretized by the forward Euler time scheme as,
As a result, the saturations are updated explicitly.
3.2 IMPES method for the discretefracture model
The semiimplicit time schemes for the pressure equations (46) and (47) in the matrix domain and fractures are formulated as,
Once the pressure equations are solved, the velocities are calculated as,
The saturations in the matrix domain and fractures are updated as,
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 , we can express the capillary pressure function given in equation (22) as a function of , i.e., . The energy parameters can be calculated by the least squares method: solve the minimizer of,
(59)where is the sample size and is the experimental value of capillary pressure measured at .
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 wettingphase residual saturations are taken as for three samples respectively, while 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 nonconvex functions with respect to the saturation , and thus the conventional convex functions [2, 4–7, 16] for the capillary pressure would not agree with the experimental data.
Fig. 1 Experimental verification of the capillary pressure formulation 
5 Numerical results
The CellCentered Finite Difference (CCFD) method is used as the spatial discretization method in numerical simulation. The detailed formulations of CCFD for the discretefracture 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 and . The absolute permeabilities in the matrix blocks and the fractures are 1 md and 10^{5} md, respectively. The quadratic relative permeabilities are applied in both fractures and matrix; that is, in equation (11). The viscosities of water and oil are cP and cP, respectively. The residual saturations are taken as and . The energy parameters are taken as bar, bar, bar in the matrix and bar, bar, 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.
Fig. 2 Example 1: distribution of fractures. 
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.
Fig. 3 Example 1: water saturation contours at different times. 
Fig. 4 Example 1: verification of thermodynamical consistency. 
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 (). Water is injected at the left end of the medium, while oil is produced at the righthand 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.
Fig. 5 Example 2: water saturation contours with free energy at different times. 
Fig. 6 Example 2: water saturation contours without free energy at different times. 
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.
Fig. 7 Example 2: free energy contours at different times. 
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.
Fig. 8 Example 3: distribution of fractures. 
Fig. 9 Example 3: water saturation contours with free energy at different times. 
Fig. 10 Example 3: water saturation contours without free energy at different times. 
Fig. 11 Example 3: energy contours at different times. 
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 energybased model has been developed for twophase flows in heterogeneous and fractured media. The model inherits the framework of the traditional model of twophase 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 energybased discrete fracture model to describe twophase 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
 Cai J., Wei W., Hu X., Liu R., Wang J. (2017) Fractal characterization of dynamic fracture network extension in porous media, Fractals 25, 1750023.
 Chen Z., Huan G., Ma Y. (2006) Computational methods for multiphase flows in porous media, Society for Industrial and Applied Mathematics, Philadelphia, PA.
 Hoteit H., Firoozabadi A. (2008) An efficient numerical model for incompressible twophase flow in fractured media, Adv. Water Resour. 31, 891–905.
 Chen J., Sun S., Wang X.P. (2014) A numerical method for a model of twophase flow in a coupled free flow and porous media system, J. Comput. Phys. 268, 1–16.
 Galusinski C., Saad M. (2008) Two compressible immiscible fluids in porous media, J. Differ. Equ. 244, 1741–1783.
 Hoteit H., Firoozabadi A. (2008) Numerical modeling of twophase flow in heterogeneous permeable media with different capillarity pressures, Adv. Water Resour. 31, 56–73.
 Hou J., Chen J., Sun S., Chen Z. (2016) Adaptive mixedhybrid and penalty discontinuous Galerkin method for twophase flow in heterogeneous media, J. Comput. Appl. Math. 307, 262–283.
 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.
 Cai J., Yu B., Mei M., Luo L. (2010) Capillary rise in a single tortuous capillary, Chin. Phys. Lett. 27, 054701.
 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.
 Sun S. (2019) Darcyscale phase equilibrium modeling with gravity and capillarity, J. Comput. Phys. 399, 108908.
 Wei W., Cai J., Xiao J., Meng Q., Xiao B., Han Q. (2018) KozenyCarman constant of porous media: Insights from fractalcapillary imbibition theory, Fuel 234, 1373–1379.
 Kou J., Sun S., Wang X. (2015) Efficient numerical methods for simulating surface tension of multicomponent mixtures with the gradient theory of fluid interfaces, Comput. Methods Appl. Mech. Eng. 292, 92–106.
 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.
 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.
 Cancès C. (2018) Energy stable numerical methods for porous media flow type problems, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 73, 78.
 Lee S., Mikelić A., Wheeler M.F., Wick T. (2018) Phasefield modeling of two phase fluid filled fractures in a poroelastic medium, Multiscale Model. Simul. 16, 1542–1580.
 Abels H., Garcke H., Grun G. (2012) Thermodynamically consistent, frame indifferent diffuse interface models for incompressible twophase flows with different densities, Math. Mod. Meth. Appl. Sci. 22, 1150013.
 Guo Z., Lin P. (2015) A thermodynamically consistent phasefield model for twophase flows with thermocapillary effects, J. Fluid Mech. 766, 226–271.
 Shen J., Yang X. (2014) Decoupled energy stable schemes for phasefield models of twophase complex fluids, SIAM, J. Sci. Comput. 36, B122–B145.
 Zhu G., Kou J., Yao B., Wu Y.S., Yao J., Sun S. (2019) Thermodynamically consistent modelling of twophase flows with moving contact line and soluble surfactants, J. Fluid Mech. 879, 327–359.
 Zhu G., Li A. (2020) Interfacial dynamics with soluble surfactants: A phasefield twophase flow model with variable densities, Adv. GeoEnergy Res. 4, 86–98.
 Chen C.Y., Yan P.Y. (2015) A diffuse interface approach to injectiondriven flow of different miscibility in heterogeneous porous media, Phys. Fluids 27, 083101.
 Cogswell D.A., Szulczewski M.L. (2017) Simulation of incompressible twophase flow in porous media with large timesteps, J. Comput. Phys. 345, 856–865.
 CuetoFelgueroso L., Juanes R. (2014) A phasefield model of twophase HeleShaw flow, J. Fluid Mech. 758, 522–552.
 Dede L., Garcke H., Lam K.F. (2018) A HeleShawCahnHilliard model for incompressible twophase flows with different densities, J. Math. Fluid Mech. 20, 531–567.
 Feng X.B., Wise S. (2012) Analysis of a DarcyCahnHilliard diffuse interface model for the HeleShaw flow and its fully discrete finite element approximation, SIAM J. Numer. Anal. 50, 1320–1343.
 Ngamsaad W., Yojina J., Triampo W. (2010) Theoretical studies of phaseseparation kinetics in a Brinkman porous medium, J. Phys. A: Math. Theor. 43, 202001.
 Sabooniha E., Rokhforouz M.R., Ayatollahi S. (2019) Porescale investigation of selective plugging mechanism in immiscible twophase flow using phasefield method, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 74, 78.
 Wise S.M. (2010) Unconditionally stable finite difference, nonlinear multigrid simulation of the CahnHilliardHeleShaw system of equations, J. Sci. Comput. 44, 38–68.
 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.
 Kou J., Sun S., Wang X. (2020) A novel energy factorization approach for the diffuseinterface model with PengRobinson equation of state, SIAM J. Sci. Comput. 42, B30–B56.
 Ghorayeb K., Firoozabadi A. (2000) Numerical study of natural convection and diffusion in fractured porous media, SPE J. 5, 12–20.
 Kazemi H., Merrill L.S. (1979) Numerical simulation of water imbibition in fractured cores, Old SPE J. 19, 175–182.
 Lee S., Jensen C., Lough M. (2000) Efficient finitedifference model for flow in a reservoir with multiple lengthscale fractures, SPE J. 5, 268–275.
 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.
 Pruess K. (1985) A practical method for modeling fluid and heat flow in fractured porous media, Old SPE J. 25, 14–26.
 Sarkar S., Toksoz M.N., Burns D.R. (2002) Fluid flow simulation in fractured reservoirs, in: Report, Annual Consortium Meeting.
 Thomas L., Dixon T., Pierson R. (1983) Fractured reservoir simulation, Old SPE J. 23, 42–54.
 Warren J.E., Root P.J. (1963) The behavior of naturally fractured reservoirs, Old SPE J. 3, 245–255.
 Amaziane B., Panfilov M., Pankratov L. (2016) Homogenized model of twophase flow with local nonequilibrium in double porosity media, Adv. Math. Phys. 2016, 3058710.
 Chen H., Kou J., Sun S., Zhang T. (2019) Fully massconservative IMPES schemes for incompressible twophase flow in porous media, Comput. Methods Appl. Mech. Eng. 350, 641–663.
 Chen Z., Huan G., Li B. (2004) An improved IMPES method for twophase flow in porous media, Transp. Porous Med. 54, 361–376.
 Coats K.H. (2003) IMPES stability: Selection of stable timesteps, SPE J. 8, 181–187.
 Fagin R.G., Stewart C.H. Jr. (1966) A new approach to the twodimensional multiphase reservoir simulator, Old SPE J. 6, 175–182.
 Kou J., Sun S. (2010) A new treatment of capillarity to improve the stability of IMPES twophase flow formulation, Comput. Fluids 39, 1293–1931.
 Kou J., Sun S. (2013) Convergence of discontinuous Galerkin methods for incompressible twophase flow in heterogeneous media, SIAM J. Numer. Anal. 51, 3280–3306.
 Kou J., Sun S. (2014) Upwind discontinuous Galerkin methods with conservation of mass of both phases for incompressible twophase flow in porous media, Numer. Methods Partial Differ. Equ. 30, 1674–1699.
 Sheldon J.W., Zondek B., Cardwell W.T. (1959) Onedimensional, incompressible, noncapillary, twophase fluid flow in a porous medium, Trans. SPE AIME 216, 290–296.
 Yang H., Sun S., Yang C. (2017) Nonlinearly preconditioned semismooth Newton methods for variational inequality solution of twophase flow in porous media, J. Comput. Phys. 332, 1–20.
 Young L., Stephenson R. (1983) A generalized compositional approach for reservoir simulation, Old SPE J. 23, 727–742.
 Cahn J.W., Allen S.M. (1977) A microscopic theory for domain wall motion and its experimental varification in feal alloy domain growth kinetics, J. Phys. Colloq. C7, C7–51.
 Chen W., Wang C., Wang X., Wise S.M. (2019) Positivitypreserving, energy stable numerical schemes for the CahnHilliard equation with logarithmic potential, J. Comput. Phys.: X 3, 100031.
 Wang X., Kou J., Cai J. (2020) Stabilized energy factorization approach for AllenCahn equation with logarithmic FloryHuggins potential, J. Sci. Comput. 82, 25.
 Kou J., Sun S. (2018) Thermodynamically consistent modeling and simulation of multicomponent twophase flow with partial miscibility, Comput. Methods Appl. Mech. Eng. 331, 623–649.
 Peng D., Robinson D.B. (1976) A new twoconstant equation of state, Ind. Eng. Chem. Fund. 15, 59–64.
 Kou J., Sun S. (2017) Efficient energystable dynamic modeling of compositional grading, Int. J. Numer. Anal. Mod. 14, 218–242.
 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).
 Kou J., Sun S., Yu B. (2011) Multiscale timesplitting strategy for multiscale multiphysics processes of twophase flow in fractured media, J. Appl. Math. 2011, 861905.
All Figures
Fig. 1 Experimental verification of the capillary pressure formulation 

In the text 
Fig. 2 Example 1: distribution of fractures. 

In the text 
Fig. 3 Example 1: water saturation contours at different times. 

In the text 
Fig. 4 Example 1: verification of thermodynamical consistency. 

In the text 
Fig. 5 Example 2: water saturation contours with free energy at different times. 

In the text 
Fig. 6 Example 2: water saturation contours without free energy at different times. 

In the text 
Fig. 7 Example 2: free energy contours at different times. 

In the text 
Fig. 8 Example 3: distribution of fractures. 

In the text 
Fig. 9 Example 3: water saturation contours with free energy at different times. 

In the text 
Fig. 10 Example 3: water saturation contours without free energy at different times. 

In the text 
Fig. 11 Example 3: energy contours at different times. 

In the text 