Thermodynamically consistent modeling of two-phase incompressible flows in heterogeneous and fractured media

. Numerical modeling of two-phase ﬂ ows in heterogeneous and fractured media is of great interest in petroleum reservoir engineering. The classical model for two-phase ﬂ ows in porous media is not completely thermodynamically consistent since the energy reconstructed from the capillary pressure does not involve the ideal ﬂ uid 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 de ﬁ ned. Consequently, the classical phase- ﬁ eld models can not be applied due to the use of diffuse interfaces. In this paper, we propose a new thermodynamically consistent energy-based model for two-phase ﬂ ows in heterogeneous and fractured media, which is free of the gradient energy. Meanwhile, the model inherits the key features of the traditional models of two-phase ﬂ ows 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 two-phase ﬂ ows in fractured media. The popularly used implicit pressure explicit saturation method is used to simulate the model. Finally, the experimental veri ﬁ cation of the model and numerical simulation results are provided.


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][2][3].
The classical two-phase flow models [2,[4][5][6][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][9][10][11][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][14][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 twophase 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 twophase flows have been extensively studied in the literature, for instance [18][19][20][21][22]. Thermodynamically consistent phasefield-based diffuse interface models and simulation of twophase flows in porous media have been reported in the literature, for instance [23][24][25][26][27][28][29][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][34][35][36][37][38][39][40], such as the single-porosity model, the dualporosity/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 twophase 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 twophase flows in fractured media.
The IMplicit Pressure Explicit Saturation (IMPES) method [42][43][44][45][46][47][48][49][50][51] is an extensively used time-stepping approach for multi-phase flow simulation.Based on multiphysics 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.

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

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 a, denoted by S a , where a = w, n, is the volume fraction of the void of a porous medium filled by phase a.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, where t is the time, u a denotes Darcy's velocity of phase a, and q a stands for the injection/production rate of phase a.
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][14][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][53][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 ra be the residual saturation of phase a and let 1 = 1ÀS rw ÀS rn .We propose the following logarithmic free energy function F as, where c a and c wn are the energy parameters that measure the capillary intensity, and S e a is the normalized saturation defined by, S e a ¼ S a À S ra 1 : 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, where g is the gravity acceleration, z is the depth, and q a is the mass density of phase a.The total energy U is the sum of the free energy and gravitational potential energy, We remark that c a and c 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 c = c w = c n = c wn , c 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: where r is the surface tension and e 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 , In heterogeneous media and fractured media, l a may be continuous in space, whereas c a and c 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 a is expressed as, The modified Darcy's law for phase a is formulated as, where p is the pressure, K represents the absolute permeability tensor of the porous medium, k ra is the relative permeability of phase a and g a is the viscosity of phase a.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 K is the symmetric positive definite tensor.The relative permeability k ra of phase a is a nonnegative function of S a , for instance [3], 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 k a ¼ kra g a and rewrite equation (10) as, 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.

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 X, we consider the simplified model, associated with the boundary conditions, where n is the outward unit normal vector to oX.
Using the constraint S w + S n = 1, we derive the divergence-free 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 k a !0 and K is the symmetric positive definite tensor, we deduce using equation ( 13) that, Thus, the proposed model is thermodynamically consistent.

Special cases
We now present several special cases of the proposed model, one of which is the classical two-phase flow model.
The a-phase pressure is denoted by p a , where a = 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 c n = c wn = 0, and as a result, l n = 0. From equation ( 22), the capillary pressure changes into p c ¼ Àl w .Darcy's velocity of the wetting phase is rewritten as, which leads to p ¼ p n as a result of, Thus, we obtain the Darcy velocity of the non-wetting phase as, From equation (7), the capillary pressure has the form p c ¼ Àc w lnðS e w Þ.So this case is the classical two-phase flow model [3,6,42,[46][47][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 c w ¼ c wn ¼ 0, and then the capillary pressure becomes p c ¼ l n .Darcy's velocity of the non-wetting phase is rewritten as, which results in p ¼ p w .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 c w ¼ c n ¼ c and c wn ¼ 0, where c is a constant.Moreover, S rw ¼ S rn ¼ 0. In this case, we have l a ¼ c lnðS a Þ and then we deduce, Taking into account equations ( 21), ( 22) and ( 28), we deduce the formulations of Darcy's velocities as, which lead to the weighted pressure as, as a result of, 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.

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, and the general capillary potential, Then we rewrite the velocity of phase a as, 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, where k t ¼ k w þ k n is the total mobility.Furthermore, we define, As a result, u t is reformulated as, We define the fractional flow function n w ¼ k w =k t .The wetting-phase velocity is rewritten as, The model is reformulated as the following governing equations, where q t ¼ q w þ q n .The boundary oX of the spatial domain X is decomposed into the Dirichlet part C D and Neumann part C N , where oX ¼ C D S C N .The boundary conditions associated to equations ( 40) and ( 41) are given by, where n is the outward unit normal vector to oX.The wetting-phase saturation on the boundary oX is subject to, The initial saturation is provided by,

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, U w;m denotes the form of U w 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, where the subscript f denotes the fracture and Q t;f stands for the total mass transfer across the matrix-fracture boundaries.The velocities u a in the matrix and fractures are expressed as, Similarly, we can formulate the saturation equation in the matrix as, and the saturation equation in the fractures as, where Q w;f represents the mass transfer of wetting phase across the matrix-fracture boundaries.

Numerical method
The IMplicit Pressure Explicit Saturation (IMPES) method [2,[42][43][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 < . . .< t M ¼ T .We denote the time step length by Át i ¼ t iþ1 À t i .The value of a function v at the time t i is denoted by v i .

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 ra are calculated from the saturations at the previous time step, and thus the variables k w , k n , and k t in the pressure equation are treated explicitly.The chemical potentials are computed explicitly as well.The pressure equation is semi-implicitly discretized as, which is a linear, symmetric positive definite equation of pressure and can be easy to be solved.Once the pressure is calculated, the velocity u iþ1 a 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.

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, Once the pressure equations are solved, the velocities are calculated as, The saturations in the matrix domain and fractures are updated as,

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 ðS w Þ.The energy parameters can be calculated by the least squares method: solve the minimizer of, where m is the sample size and 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 rw ¼ 0:098; 0:039; 0:03 for three samples respectively, while S 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][5][6][7]16] for the capillary pressure would not agree with the experimental data.

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 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 / m ¼ 0:2 and / f ¼ 1.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, l w ¼ l n ¼ 2 in equation (11).The viscosities of water and oil are g w ¼ 1 cP and g n ¼ 0:5 cP, respectively.The residual saturations are taken as S rw ¼ 0:03 and S rn ¼ 0. The energy parameters are taken as c m w ¼ 0:1166 bar, c m n ¼ 0:0108 bar, c m nw ¼ 0:0742 bar in the matrix and c f w ¼ 1:166e À 3 bar, c f n ¼ 1:08e À 4 bar, c f nw ¼ 7:42e À 4 bar in the fractures.The effect of gravity is neglected in the horizontal media.

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.

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

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.

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.

Fig. 1 .
Fig. 1.Experimental verification of the capillary pressure formulation

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

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