Open access peer-reviewed chapter

Porous Flow with Diffuse Interfaces

Written By

Paul Papatzacos

Submitted: September 7th, 2020 Reviewed: December 12th, 2020 Published: January 19th, 2021

DOI: 10.5772/intechopen.95474

Chapter metrics overview

266 Chapter Downloads

View Full Metrics


This chapter presents a model developed by the author, in publications dated from 2002 to 2016, on flow in porous media assuming diffuse interfaces. It contains five sections. Section 1 is an Introduction, tracing the origin of the diffuse interface formalism. Section 1 also presents the traditional compositional model, pointing out its emphasis on phases and questioning the concept of relative permeabilities. Section 2 presents the mass, momentum, and energy balance equations, for a multicomponent continuous fluid, in their most general form, at the pore level. The existence of constitutive equations with phase-inducing terms is mentioned, but the equations are not introduced at this level, and phases are not an explicit concern. Section 3 is about the averaging of the pore level equations inside a region containing many pores. There is no explicit mention of phases and therefore not of relative permeabilities. Section 4 is the technical basis from which the constitutive equations of the model arise, and it is shown that many models can exist. Section 5 introduces constitutive equations and presents a minimal model for multicomponent, multiphase, and thermal flow in neutrally wetting porous media, i.e., a model with a minimal amount of phenomenological parameters.


  • flow in porous media
  • Marle averaging
  • diffuse interface
  • multiphase flow
  • phase segregation
  • relative permeabilities

1. Introduction

The model presented in this chapter was developed in Refs. [1, 2, 3]. Ref. [1] covers the one-component two-phase case, Ref. [2] is a generalisation to an arbitrary number of components and three phases (two liquids, one gas), Ref. [3] generalises to variable temperature. All three show applications to neutrally wetting media. The way wetting can be accounted for is discussed in general terms in Ref. [1]. A partial practical implementation, valid for incomplete wetting, is suggested in Ref. [4] but a fully satisfactory solution that accounts for total wetting (or capillary condensation in adsorption terms) is still a matter for further research.

The diffuse interface theory was initiated in 1893 in a paper by van der Waals (see the translation by Rowlinson [5]) where he proposes to replace the old assumption of a surface between phases by the assumption of a continuous transition inside a thin interphase region, where certain quantities, notably the density, vary continuously. The core of his theory consists of a Helmholtz function modified by the addition of a term proportional to the squared gradient of the density, thus accounting for the energy stored in the region. In the model presented in the following pages, the van der Waals theory is introduced at the upscaled level. See Section 4. Generalising the van der Waals expression to ν chemical components [2], one obtains


where Rα is the upscaled density of fluid component α, F is the Helmholtz function per unit volume, the superscript F referring to the fluid, while bF means “bulk fluid” and refers to fluid regions that are far away from interphase layers. The Λα are constants (Assumption A7, Appendix 1). For later reference FbF and its differential are


where PbF is the pressure in the bulk fluid, SbF the entropy per unit volume, and CbFα the chemical potential of component α, divided by its molar mass (the Cα/Mα of Refs. [2, 3]).

The van der Waals paper was followed, in 1901, by a paper by Korteweg [6] about the equations of motion of a fluid with large but continuous density changes, where Korteweg showed that, for such a fluid, the usual scalar pressure must be replaced by a symmetric second order tensor.

The van der Waals-Korteweg papers were apparently forgotten, then rediscovered in the nineteen seventies, the diffuse-interface method being introduced as a novel way to solve fluid mechanics problems in two-phase flow. For a review, see [7]. While van der Waals and Korteweg assumed that the gradients in Eq. (1) are small, modern advances have shown that this limitation can be lifted [7, 8].

Ref. [8] is, to the author’s knowledge, the first formulation of flow in porous media with diffuse interfaces.1 The purpose with such a novel formulation is to avoid some of the known weaknesses of the traditional compositional models treating multiple phase problems. The mathematical core of the compositional models used in reservoir engineering, consists of equations expressing mass balance for each chemical component. The distribution of the phases in the reservoir is essential to the formulation, and is determined at the beginning of each time step. A component, present in a phase, is transported with the Darcy velocity of the phase. Each phase-dependent Darcy velocity is written with a permeability that is modified by a multiplicative factor, the relative permeability. The above description emphasises two of the weaknesses of the traditional models. The first is the assumption that well-defined phases exist at all times. The second is the use of relative permeabilities, a concept that is “seriously questioned”, as expressed by Adler and Brenner in a 1988 paper [9]. Concerning the understanding of relative permeabilities in the framework of the model presented here, see the comments following Eq. (100) below.

The mathematical core of the model (actually a family of models) to be presented consists of mass balance equations, one per chemical component, of a momentum balance equation, and of an entropy balance equation. The thermodynamical description of the fluid mixture involved is part of the core. The central purpose is to calculate the component densities, and other characteristic quantities such as fluid velocity and fluid temperature, as functions of space and time. If the approximation of constant temperature is valid, only the mass and momentum balance equations are necessary.

Phases (and thereby relative permeabilities) do not take part in the formulation. They result from the solutions of the model equations, and are detected by rapid variations of densities, and by regions of approximately uniform densities. They can be shown to exist in static equilibrium or steady state dynamical situations [1, 2, 3].

A minimal model is presented in Section 5, consisting of a minimal amount of parameters.

A note on the appendices: Some concepts are grouped in appendices for easy reference. Appendix 1, for example, lists all the assumptions the model is built on.

A note on wetting: The problem of accounting for the wetting properties of the pore surface remains to be solved. There are two approaches to the problem: through boundary conditions to the Navier–Stokes equations at pore level, or through the theory of adsorption at the upscaled Darcy level. The first approach has been used in the publications considered here, and it is explained in Ref. [10] that the diffuse interface theory presented in Ref. [1] is consistent for neutral wetting, i.e., for pore level wetting angles around 90: see Appendix 1, Assumption A0. Refs. [2, 3] assume this limitation.

A note on notation: Right-handed Cartesian coordinates x1x2x3 are assumed, the plane x1x2 being horizontal, and axis x3 pointing upwards. Any vector A has components Ak where k is 1, 2, or 3. In addition, the notation t=/t and k=/xk is used. The summation convention applies to latin indexes i, j, k, l: an index that is repeated in a term (as in AjBj, or Ckk) indicates summation (j=13AjBj, or k=13Ckk). Note that an index that is repeated, but not in the same term (as in Ai=Bi), means that the expression is valid for all values of the index in the set 1,2,3. Symbols are otherwise defined when introduced.


2. Pore level equations

The fluid is a mixture of ν chemical components, it is continuous, and if phases exist, there are interphase regions where quantities vary continuously, possibly rapidly. With assumptions A1 to A3 (Appendix 1) the balance equations for mass, momentum, and energy, are written below in a most general manner:


For the theoretical basis of these equations see chapter 11 of the book by Hirschfelder et al [11]. Greek superscripts indicate the species so that ρα is the mass density of component α, while ρ is the total density:


Further, v is the local velocity, defined as the total momentum divided by the total mass, both per volume; ια is the non-convective mass current of component α, with the property


J is the non-convective energy current, and tij is the stress tensor; ε=u+12ρv2 is the total energy, u is the internal energy, both per unit volume, and f is an external force per unit mass of fluid, thus being species-independent. In the case of gravity


It is convenient, for later use, to introduce the gravitational potential


The symmetry tij=tji is assumed below. It is reminded that it expresses angular momentum balance for a fluid where no other torques exist than the one due to the external force f, and the one due the surface stress nt acting on every surface element where the normal vector is n.

Within the van der Waals theory, one expects that ιkα, Jk, and tki, contain terms whose magnitudes are important in the interphase regions, but are otherwise negligible. No constitutive equations are introduced at the pore level, but it is mentioned for later reference, that for a simple one-phase fluid,


where p is the pressure and θij is the viscous stress tensor (symmetric and linear in the gradients of the vi, with coefficients of shear and bulk viscosity).

The upscaling, i.e., the averaging over many pores, is done in the next section by the method due to Marle [12]. This method assumes that the physical quantities that appear in the balance equations above are treated as distributions [12, 13], the underlying reason being that such quantities are discontinuous, and that one needs to average their partial derivatives. Taking ρ as an example, it is (i) a fluid density in a pore, (ii) a rock density in the rock matrix, (iii) undefined on the pore surface, and one needs to average tρ and kρ.

The physical quantities appearing in Eqs. (3) to (5) are listed in the first line of Table 1. Their values in the pores, or in the rock, are denoted with a superscript F, or S, as shown in the second and third lines of the table. In the third line, a missing entry indicates non-existence, the first two 0-values indicate no material transport in the solid, the third zero value follows from Assumption A4 (Appendix 1): as shown by Marle [12], this assumption implies that the momentum balance equations in the solid and on the pore surface have the form 0=0. It is important to keep in mind, especially when averaging, that any quantity with an F superscript is equal to 0 on the pore surface (denoted by Σ) and in the solid; likewise, any quantity with superscript S is equal to 0 on Σ and in the fluid. The quantities in the first line of the table are not defined on Σ, but limit values are, as in

In fluidρFαρFvkFιkFαtkifkuFjkF
In solidρS000uSjkS

Table 1.

First line, left of the vertical: quantities appearing in Eqs. (3) and (4); right of the vertical: quantities only appearing in Eq. (5). Second and third lines: notation when specialising to the fluid and the solid. Concerning the missing entries and the three 0-values, see beginning of paragraph containing Eq. (11).


expressing that the fluid velocity vanishes on Σ. Expressions of the type XΣ for some X, used in this section and in the next, define the limit of the quantity X at Σ, along the line carrying the normal to Σ pointing towards the fluid. The existence of the normal implies some idealisation of the pore surface.

Note finally that, according to Eq. (6),


As a preliminary to upscaling, Eqs. (3)(5) are now explicitly written in terms of distributions, using the notation of Appendix 2 where A is replaced by Fα or F, as appropriate, while B is replaced by S. As an example, ραxt, is a distribution depending on space x and time t. It is continuous in time but not in space, being equal to ρFαxt if x is inside a pore and to ρSxt if x is in the rock. It is not defined when x is on Σ but we assume that ρFαΣ and ρSΣ exist.

The generalised mass balance equation is obtained from Eq. (3), using Eqs. (112), (113), and (11):


From this equation one now gets, using Table 1, the equations that are separately valid inside the pores, on Σ, and in the solid:

ιFαnΣδΣ=0,      α=1ν,onΣE15

Turning to the generalised momentum balance equation, one must account for Assumption A4 (Appendix 1), implying that the momentum balance equations in the solid and on the pore surface have the form 0=0. Inside the pores, one simply re-writes Eq. (4), with the superscript F on the density and velocity:


The generalised energy balance equation is obtained from Eq. (5), using Eqs. (112), (113), and (11):


Using Table 1, one then obtains:


3. Averaging

The Marle averaging process [12] is followed in all essentials, except in the assumption that well-defined phases exist, separated by interphase surfaces. The averaging volume is a sphere of radius r, large when compared to a pore radius, small when compared to a linear dimension of the reservoir. A C function mx is introduced, somewhat flat around x=0 and equal to zero for xr, normalised so that its integral over all space is equal to 1. (See Eq. (17) in Ref. [12] for an example of such a function.) Given any function of space and time, fxt, its average Fxt is obtained by the convolution


where the integration is over all of space. The convolution ensures that F is C [12].

The averaged balance equations are differential equations in the averaged quantities (averaged densities, velocities, ...). These equations are established by a three-step process. Step 1: the generalised equations for mass, momentum, and energy balance are each in turn convoluted with m. Step 2: the following rules [12] are applied, allowing to take the differential operators out of the averaging convolutions:


where fZfZxt and Z is F, Fα, or S; εZ is 1 if Z is F or Fα, −1 if Z is S. Eq. (11) is also applied at step 2. Step 3: the remaining convolutions are used to define averaged quantities, where the following constraints should be obeyed: (i) except for the definition of porosity, the way to define the averages is suggested by the equations obtained after completion of step 2; (ii) the averaged equations have essentially the same forms as Eqs. (3) to (5).

Differential equations in the averaged quantities result from the three steps. It is shown below that the mass and momentum balance equations should be treated together, and that the energy equation can be treated as an addition.

The averaging of the mass and momentum balance equations follows. Steps 1 and 2 are applied to Eqs. (14) to (17), and lead to:


Note that the last term on the right-hand side of the last equation originally is ρFfim, but fi can be taken out of the convolution integral because of Assumption A5 (Appendix 2). The definitions of averaged (or upscaled) quantities now follow (step 3). Porosity Φ is defined first, then follow definitions suggested by the convolution operations in the equations above.

Porosity Let a function χx be 1 when x is in F, and let it be 0 otherwise.



Note that, according to definition (22), Φ can depend on x.

Species adsorption The left-hand side of Eq. (26) defines the amount, KΣα, of component α adsorbed at the pore surface:


Density of solidEq. (27) suggests defining the solid density RS by:


Density of componentα The first term on the left-hand side of Eq. (25) suggests defining the density Rα of fluid component α by


Note that Eq. (12) implies that the averaged total fluid density, R, is


Fluid velocity The first term on the lef-hand side of Eq. (28) has the convolution of a product of two term, where the average of one of them is known from Eqs. (32) and (33). The averaged fluid velocity, denoted Vi (without the F superscript since there is no velocity in S or on Σ) is then defined by


Diffusive mass current of componentα The second convolution on the left-hand side of Eq. (25) is used to define the upscaled diffusive current in the fluid, denoted Ikα (without the F superscript since there are no diffusive currents in S or on Σ). It is not defined as the average of ιkFα because both pore level effects of convection and diffusion contribute to it [12]. Keeping in mind the constraint that the averaged equations should have the same form as the original ones, Ikα is defined by


Note that, summing this equation over α from 1 to ν, and using previous results, one gets


Stress tensor The second convolution on the left-hand side of Eq. (28) suggests defining Tki, the upscaled version of tki, by


Note that tki=tik implies Tki=Tik.

Frictional force per unit volume The first convolution on the right-hand side of Eq. (28) suggests defining the frictional force per unit volume FiF by


The upscaled mass and momentum balance equations now follow, in the following order: mass balance for the solid, mass balance at the pore surface, mass balance for the fluid in the pores, and momentum balance for the fluid in the pores:


The first equation states that porosity and solid density do not vary with time (consistently with Assumption 4 (Appendix 2) but can vary in space. The second equation states that adsorption is negligibly small, for any component, consistent with Assumption A0 (Appendix 2). The remaining two equations determine the ν component densities and the three velocity components. This means that the components of the diffusive mass current, of the stress tensor, and of the frictional force must be provided. (Constitutive equations).

The averaging of the energy balance equations now follows. Steps 1 and 2 are applied to Eqs. (19) to (21), giving:


As in the case of Eq. (28), Assumption A4 (Appendix 2) has been used to take fk out of the convolution on the right-hand side of Eq: (43). Six definitions are introduced below, built on the definitions that were introduced in connection with the averaging of the mass and momentum balance equations. An underscore indicates the defined quantity.

Internal energy per unit volume of solid


Internal energy current of solid


Solid to fluid energy transfer


Fluid to solid energy transfer


Internal energy per unit volume of fluid


Internal energy current of fluid


Using these definitions in Eqs. (43) to (45) one obtains:


Eq. (52) contains a redundancy in the form of a balance equation for kinetic energy. This equation can be obtained directly by multiplying both sides of Eq. (42) with Vi, summing over i, and using Eq. (41). Subtracting the equation thus obtained from Eq. (52), one gets


4. Basis for constitutive equations

In practical application, the upscaled balance equations will be used to calculate, primarily, the densities Rα, the velocity components Vi and the temperature T. Assuming that the approximation of constant T is valid, one needs only focus on the mass and momentum balance equations, (39) to (42). One sees in this case that one needs expressions for the Ikα, the Tij, and the Fi in terms of T and the Rα. If the approximation of constant temperature is not valid, one must first distinguish between the temperatures of the solid and the fluid, and one needs an equation for the energy transfer in case of a temperature difference. One introduces the simplifying Assumption A6 (see Ref. [12] and Appendix 1), from which it follows that just one additional equation is needed for calculating Txt. Such an equation usually describes the evolution of either total energy or total entropy. In either case, expressions are needed for the currents JkS and JkF. (Expressions for the energy transfers QSF and QFS are unnecessary since they cancel when taking the sum of the solid and fluid energies). Most of what is needed is obtained in Section 5 by applying the theory of irreversible processes, starting from the evolution equation for entropy, although it seems that some preliminary work is unavoidable to directly obtain an expression for the pressure tensor Pij that replaces the usual scalar pressure.

The derivation of the pressure tensor is given below, followed by the derivation of the evolution equation for the total entropy. It is essential to useexpression (1)in both derivations.2

4.1 The pressure tensor

One considers the upscaled fluid, consisting of a mixture of ν components in a container with surface ∂Ω and volume Ω, inside a large bath at uniform and constant temperature T. One looks for conditions of equilibrium in the presence of gravity. The fluid has the Helmholtz free energy density given by Eq. (1), and it is assumed that the bounding surface is neutrally wetting so that there is no energy stored on ∂Ω. The total energy stored in the fluid is


where W is the gravitational potential (see Eqs. (8) and (9)).

One now looks for the conditions the Rα satisfy when F is at its minimum, given that the total mass, ΩR, is constant. This is minimisation with constraint, a standard problem in variational calculus. It is easily found that


where κ is a Lagrange multiplier, possibly a function of T but not of the Rα. See Ref. [14] for the technique and the theorems involved: essentially, the expression on the left-hand side of Eq. (57) must be continuous. Note that Eq. (58) expresses the non-wetting property of the outer boundary: the density neither increases nor decreases along the normal.

The next step consists in multiplying Eq. (57) with iRα and summing over α. Each term of the resulting equation can then be re-written as a gradient (first term), or as a sum of a gradient and the component of a force (sum of second and third terms), or as a sum of a gradient and a divergence (fourth term). The result is the following expression:




The first equation obviously generalises the classical P=f, and it is natural to generalise Eq. (10) by setting


Tij and Pij are symmetric, so that Θij=Θji. Eqs. (60) and (61) are used in the long calculation that lead to the evolution equation for the entropy. They are also used to get expressions (99) and (100) for the modified Darcy equation, where the presence of the Λ‘s is exclusively due to their presence in Eq. (60).

4.2 The entropy evolution equation

Eqs. (1) and (2), together with Eqs. (41), (54), (55), (60), and (61), are needed to get the entropy equation for the fluid.

Taking the differential of Eq. (1) and using Eq. (2), one obtains, keeping in mind the assumption that the Λα are constants,


This equation leads immediately to two conclusions. There is no additional entropy, and no additional chemical potentials due to large density gradients, since: (i) SF, defined as FF/T, is equal to SbF; (ii) CFα, defined as FF/Rα, is equal to CbFα. The differential of FF can now be re-written with a simpler notation:


To obtain dSF in terms of dUF and dFF one differentiates UF=TSF+FF, obtaining


This expression is now used to construct an evolution equation for the total entropy (including fluid and solid). Such an equation must be of the form


where SS is the entropy of the rock, Pk is a diffusive current, and Q is a source term.

Eq. (64) directly gives the two equations that follow, by replacing d with t (first equation), then with Vii (second equation):


Since it is ΦSF that is needed to get Eq. (65), one needs to multiply the two equations above with Φ and then commute Φ with t and i. The same commutations are necessary on the right-hand sides of the above equations, so as to get ΦUF and thus allow the use of Eq. (55). Assumption A8 (Appendix 2) has been added to Assumption A4 so as to avoid the proliferation of iΦ-terms in the minimal model to be presented. With this simplification in place, and remembering to replace Tij in Eq. (55) with the expression obtained from Eqs. (60) and (61), one obtains, after elementary but somewhat long calculations:


Keeping in mind that one is looking for an entropy equation of the form of Eq. (65), and in anticipation of using the methods of irreversible processes, one now considers each term, or group of terms, on the right-hand side above and writes it either as a divergence, or as the sum of a divergence and a scalar product. Most terms on the first line already have the required form. The general term in the first sum on the second line is easily transformed as required, the general term in the second sum also, although with somewhat more work. As to the third line, it is easily seen to be a sum of divergences. One then gets:

TtΦSF+iΦSFVi= kJkF+Φα=1νAαIkα+Φα=1νΛαkRαRαiVi+iIiα+QSFFkFVk+ΘjijViΦα=1νIiαiAα,E69

where Aα=CFαΛα2Rα. Note that Aα occurs twice on the right-hand side of Eq. (69), once as Aα, and once as iAα, both times multiplying Ikα in a sum over α. Eq. (36) then implies that one can modify the above expression of Aα by the addition of any expression that does not depend on α. To conform with the notation of Ref. [2], and especially Ref. [3], one then sets


as the expression to substitute on the right-hand side of Eq. (69).

The entropy equation for the solid is much easier to obtain because of Assumption A4 (Appendix 2). Indeed, Eq. (64) is replaced by


which implies


and, using (54):


Taking the sum of Eqs. (69) and (73) one sees that QSF cancels, and that the right-hand side consists of a sum of scalar products and of a divergence, say kKk. A last step remains because the left-hand side has a multiplicative factor T (see Eqs. (69) and (73)). Dividing both sides by T does not affect the scalar products but it modifies the divergence, at least if one assumes that T is variable, to kKk/T. This is easily remedied by writing


One finally obtains




QFS and QSF do not appear in the entropy equation, and are thus not determinable inside the model. Similarly, the currents JkF and JkS only appear as summed, so that they are not determined individually inside the model. (See the next section.)


5. Constitutive equations and the minimal model

The source term in the entropy equation plays a central role in what follows. It has been written as a sum of scalar products. Each term of this sum is the scalar product of a force (explicit or generalised) and a current. FkF is an explicit force, while the gradient of some quantity (temperature, velocity component, ...) is a generalised force. The theory of irreversible processes states that linear relations exist (at least for processes not far from equilibrium), between forces and currents. The coefficients, called phenomenological coefficients, are parameters whose signs must be such that the source term cannot be negative, to ensure against a decrease of the entropy when the system is isolated. The linear relations just mentioned are constitutive equations, and it is implied that the phenomenological coefficients must be provided as input. For various reasons, some coefficients can be put equal to zero, one often recurring reason being the belief that they are negligible in the physical situation considered. Thus there is a family of models, each member being characterised by its set of non-zero phenomenological coefficients.

The previously mentioned minimal model (see the text after Eq. (67)), is the model that contains the least possible number of non-zero phenomenological coefficients. “Least possible” means that one must obey the constraints that exist for these coefficients (Onsager symmetry, isotropy, sign) and arbitrarily setting some of them equal to zero is not always possible. A certain amount of trial and error is also required to avoid unduly reducing the model’s predictive power.

The usual vector notation is used for tensors of order 1, i.e., a bold faced letter is used when the subscript can be suppressed. Vector forces are here denoted X (with components Xi), and vector currents are denoted Y (with components Yi), with a superscript to discriminate between the different currents and forces in an alphabetic order aa shown below. (Y has been chosen instead of the traditional J so as to avoid confusion with the currents related to the energy equations.) The usual tensor conventions are assumed for latin subscripts (see “A note on notation” in the Introduction). Concerning the second order tensors Θji and jVi, on the right-hand side of Eq. (77), one sets ΔjijVi and, referring to Appendix 3, especially to Eq. (115), one writes


where XA and XijB are found by replacing Z by Δ in Eqs. (115). Note that YijB, being symmetric, has no antisymmetric part.

Referring now to the vectors on the right-hand side of Eq. (77), one introduces the following notation:


Using Eq. (116), one easily finds that


Note that Q is a proper scalar since the source of entropy is independent of the coordinate system, and does not change sign when the coordinate system changes handedness. The right-hand side of Eq. (81) is then a sum of scalar products of proper tensors, the only pseudo vector, X˜k, having dropped out.

Referring to the first sentence of this section, one writes the currents as linear combinations of the forces:


In these expressions, the L are the phenomenological coefficients. They are independent of the generalised forces, but they can depend on the temperature and the component densities. They are tensors, their orders being equal to the number of subscripts. They obey the Onsager relations [11, 15]: any coefficient with n subscripts, denoted say, LnMN, obeys


where the n subscripts are the same but not necessarily in the same order. Details concerning subscripts are not needed in what follows.

5.1 The viscosity tensor Θij

With hindsight, one knows that the upscaled viscosity tensor is not required, so that one is justified in setting equal to zero all the L‘s in the first two lines of the system of Eqs. (82), and also setting to zero all the L‘s related to the zeroed ones by the Onsager symmetry:


Then YA=YijB=0, so that


and the source term of the entropy equation reduces to


5.2 Vector currents and forces

System (82) now reduces to linear relations between vectors, the coefficients being second order tensors:


where Onsager symmetry is accounted for. Remembering that Yiα=Ikα, Eq. (36) leads to three constraints on the L-coefficients of Eq. (89):


The linear combinations above show the possibilities of constructing models where interactions between thermal conduction, fluid flow, and mass diffusion are quantified by choosing the L coefficients. Also, thermal conduction and permeability can be modelled by second order tensors, through tensors LikCC and LikDD.

However, if one limits oneself to the minimal model where: (i) the diffusive entropy current, Ki/T, is only due to the kT force, (ii) the FkF frictional force is only due to fluid velocity, Vi, (iii) each non-convective mass current, Iiα, is only due to the gradients of the Aβ, then one requires


not violating the constraints in Eqs. (90).

Note that it is not possible to simplify the model to the extent that all cross-couplings are eliminated, since that would imply that Likαβ are zero except when the superscripts are equal: such a matrix would not obey the third constraint in (90). One could of course set Likαβ=0 for all α and β, but that would eliminate all the non-convective mass currents from the model, and probably make it useless.

In the minimal model one can add a fourth requirement to the three above: the upscaled fluid and the upscaled medium are isotropic. Then further simplifications result since the remaining second order tensors, LikCC, LikDD, and Likαβ, that are properties of the solid and the fluid, must be invariant under rotations of the coordinate axes. Such tensors are called isotropic and it can be shown that an isotropic second order tensor is proportional to the Kronecker delta (see Ref. [15]). Isotropy thus introduces the following restrictions:


and one gets


The source term of the entropy equation is now


Keeping in mind that the αβ-matrix is symmetric, and that the sum of all elements on the same line or the same column is zero, it can be shown that


According to the expression for the source of entropy above, the remaining phenomenological coefficients must then satisfy


These coefficients are determined below.

Using the notation of expressions (80) in the first of expressions (93) one gets Ki proportional to iT. Now Ki/T is an entropy transport by diffusion, so that Ki is heat transport by diffusion. To recover Fourier’s law one sets


where k>0 is the thermal conductivity of the averaged solid–fluid system. See Ref. [3].

The second of expressions (93) gives FiF=Vi/DD, and it is shown below that the choice


away from interphase regions. In this expression, η is the viscosity of the pore fluid, and K is the absolute permeability. As already mentioned, the phenomenological coefficients can be functions of the component densities and of the temperature, and it is known that η is such a function (see Section 3.4 in Ref. [3]). The possibility of letting K be such a function is not used in the minimal model. To prove the implication (98), one uses Eq. (42), with the upscaled stress tensor given by Eqs. (60), (61), and (85). Using Eq. (41) one obtains:


The sum over α vanishes in the bulk fluid and one is left with the Darcy formula with an additional term, proportional to the material derivative of the velocity. One can carry out order of magnitude estimates of the three first terms in the square brackets above, in the manner of Section 4.1 of Ref. [3]. Using the numerical values given in the Appendix of the same reference, one easily finds that, if the gradient of pressure is of order 1, then the gravity term Rfi is of order 103, and the material derivative term is of order 109. Neglecting the material derivative, one obtains the Darcy formula, modified by terms that only become significant inside the interphase regions. Specialised to a one component fluid, this modified Darcy formula is:


It is shown in Ref. [4] that the added non-Darcy term can, in some well-defined flow types, produce a relative permeability when its numerical contribution is taken away as an added term, then put back as a multiplicative factor to the Darcy term. However, it is concluded in Ref. [4] that relative permeabilities cannot capture the full complexity of two-phase flow.

Given below is another version of Vi, that follows by application of the Gibbs-Duhem equation, obtained by differentiating FbF and using the expression of dFbF (Eqs. (2)):


where Aα is given by Eq. (70). Using now the notation of expressions (79) in the third of expressions (93) one gets the diffuse mass current of component α:


where it is reminded that the αβ-matrix is symmetric, its off diagonal elements are negative, and the elements of any line (or column) sum to zero. In addition, the elements can be functions of the densities. This last turns out to be extremely useful because the kCFβ behave as 1/Rβ when Rβ goes to zero, which is not acceptable when the differential equations of the model are solved numerically. In the minimal model one then sets


using the same positive number a for all elements. Each diagonal element αα is then the negative sum of the off diagonal elements of line α. For an example of such a matrix in the case ν=3, see Refs. [2, 3].

5.3 Closing details

It is easy to see that the first two restrictions displayed in Eqs. (92) do not implicate any other assumptions done in the minimal model so that the restrictions can be lifted, either singly or together, thus allowing thermal conductivity and/or permeability to be represented by a second order tensor when experiments indicate that such upgrading is required.

A non-thermal version of the model consists of ν mass balance equations, see Eqs. (41), where the Darcy velocity is given by Eq. (101) and the mass diffusion velocities by Eqs. (102) and (103), the auxiliary variables being defined by Eqs. (70). Fot the thermal version of the model one needs to re-write the entropy equation, Eq. (75), in terms of temperature: see Section 2.3 of Ref. [3].

Start and boundary conditions must be supplied for the numerical solutions of the differential equations of the model. Special attention must be taken with the boundary conditions since the equations are of the fourth degree in the space variables. See Ref. [3] for a detailed presentation.

The transport coefficients, k and η, of Eqs. (97) and (98), are needed as functions of T and of the Rα. See Ref. [3] and references given there.

The central thermodynamical function of the model is the Helmholtz function, especially in the bulk, introduced by Eq. (1). It is calculated from the equation of state of the mixture considered, which must be van der Waals or related (Redlich-Kwong, ...) so that, for temperatures less than the critical, regions of unstable fluid insure the existence of interphase regions; association terms must be included for the polar molecules of the mixture (see Ref. [16] and references given there).

It is shown in Ref. [16] that, independently of the equation of state that is chosen, the Helmholtz function contains the sum of the Helmholtz functions of the components, each considered as a gas where molecular interactions are neglected (ideal). Each ideal gas Helmholtz function contains a function of T that drops out under differentiation when T is assumed constant, but is important for the thermal model in accounting for the energy stored in the internal degrees of freedom of the molecules. The Helmholtz function of one mole of an ideal gas of component α is (see Refs. [2, 16] and references given there):


where vα is the molar volume, RG is the gas constant, e=exp1, and VαT is a function of temperature with the dimension of a volume per mole. Vα depends on the atomic masses, the principal moments of inertia of the molecule, .... It can be obtained from statistical physics, as shown in Ref. [16].

It is briefly shown below that Vα can also be obtained from an expression of fid.gas,α in terms of the heat capacity of the component. Assuming Pvα=RGT, one easily finds that, with T and vα as independent variables, the molar internal energy uα (the “id.gas” superscript is dropped for simplicity) is independent of vα, and that one can write


where cVα is the molar heat capacity of component α. The molar entropy, sα, obeying Tdsα=duα+Pdvα, one obtains


Integrating the last two displayed equations, one obtains fid.gas,α as uαTsα:


where u0α and s0α are the internal energy and the entropy at a reference state where volume, temperature, and pressure are v0α, T0, and P0=RGT0/v0α. Equating the right-hand sides of Eqs. (104) and (107), one obtains an expression for TlnVα in terms of the experimentally measurable function cVα. As stated in Ref. [16], it is in fact the derivative with respect to T of TlnVα that is needed in the differential equations of the thermal model. Assuming s0α=0, one obtains


Concerning examples of numerical solutions of the equations of the minimal model, see Ref. [2] for phase segregation, and for coning at uniform temperature; see Ref. [3] for an injection-production situation at variable temperature.



A.1 Appendix 1: Assumptions

A0: None of the chemical species completely wets the rock

A1: There are no sources or sinks

A2: There is no loss of energy by radiation

A3: There are no chemical reactions between the chemical components

A4: The solid is perfectly rigid

A5: The external force per unit volume fi is approximately constant in the averaging volume

A6: At each point, the difference between the solid and fluid temperatures is negligible

A7: The Λα are constant numbers inside the porous medium

A8: The porosity Φ is uniform in space


A.2 Appendix 2: Derivatives of distibutions – Formulas

The following is a set of formulas for the space derivatives of distributions having discontinuities across a given surface Σ. For proofs, see [12, 13].

Consider a function fxt symbolising a physical quantity, continuous in time but dicontinuous in space across a surface Σ that divides space in two regions, called A-side and B-side; f=fB on the B-side, f=fA on the A-side. Σ has a normal vector n at each point, pointing towards A:


One now defines the following regular distributions:


Then [12, 13].


where δΣ is a surface-Dirac distribution, its action δΣφ on a so-called test function φ [13] being defined as:


A.3 Appendix 3: Formulas concerning tensors

Let Zij be an arbitrary second order tensor. It can be shown (see Ref. [15]) that the following expression has general validity:


If Zij is a tensor, then Z is a scalar, Ẑ is a symmetric traceless tensor, and Z˜k is a pseudo-vector. A most useful property is as follows: if Uij and Vij are two tensors then



  1. 1. Papatzacos P: Macroscopic two-phase flow in porous media assuming the diffuse-interface model at pore level. Transport in Porous Media 2002;49:139–174. DOI: 10.1023/A:1016091821189
  2. 2. Papatzacos P: A model for multiphase and multicomponent flow in porous media, built on the diffuse-interface assumption. Transport in Porous Media 2010;82:443–462. DOI: 10.1007/s11242-009-9405-2
  3. 3. Papatzacos P: A model for multiphase, multicomponent, and thermal flow in neutrally wetting porous media, built on the diffuse-interface assumption. Journal of Petroleum Science and Engineering 2016;143:141–157. DOI: 10.1016/j.petrol.2016.02.027
  4. 4. Papatzacos P, Skjæveland SM: Diffuse-interface modelling of two-phase flow for a one-component fluid in a porous medium. Transport in Porous Media 2006;65:213–236. DOI: 10.1007/s11242-005-6081-8
  5. 5. van der Waals JD: The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density, Journal of Statistical Physics 1979;20:200–244. DOI: 10.1007/BF01011514
  6. 6. Korteweg DJ: Sur la forme que prennent les equations etc, Archives Neerlandaises des Sciences Exactes et Naturelles, 1901;II:6–20. DOI: Not known
  7. 7. Anderson DM, McFadden GB, Wheeler AA: Diffuse-interface methods in fluid mechanics. Annual Review of Fluid Mechanics 1998;30:139–165. DOI: 10.1146/annurev.fluid.30.1.139
  8. 8. Ganesan V, Brenner H: A diffuse-interface model of two-phase flow in porous media. Proceedings of the Royal Society London A 2000;456:731–803. DOI: 10.1098/rspa.2000.0537
  9. 9. Adler PM, Brenner H: Multiphase phase flow in porous media. Annual Review of Fluid Mechanics 1988;20:35–59. DOI: 10.1146/annurev.fl.20.010188.000343
  10. 10. Papatzacos P, Skjæveland SM: Relative permeability from thermodynamics. Society of Petroleum Engineers Journal 2004;9:213–236. DOI: 10.2118/87674-PA
  11. 11. Hirschfelder JO, Curtiss CF, Bird RB: Molecular theory of gases and liquids, John Wiley and Sons, Inc. New York, Chapman and Hall, Limited, London, 1954. DOI: 10.1002/pol.1955.120178311
  12. 12. Marle CM: On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media. International Journal of Engineering Science 1982;20:643–662. DOI: 10.1016/0020-7225(82)90118-5
  13. 13. Appel W: Mathematics for physics and physisists, Princeton University Press 2007. ISBN-13: 978–0–691-13102-3
  14. 14. Weinstock R: Calculus of Variations, Dover 1974. ISBN-13: 0–486–63069-2
  15. 15. Reichl LE: A Modern Course in Statistical Physics, 2nd Edition, John Wiley and Sons 1998. ISBN: 0–471–59520-9
  16. 16. Papatzacos P. The Helmholtz free energy of pure fluid substances and fluid mixtures. In: Proceedings of the 8th International Conference on Heat Transfer, Fluid Mechanics and Thermodynamics, (HEFAT2011); 11–13 July 2011; Mauritius


  • The author first became aware of this paper in 2015. There is one important difference with what is presented in the pages that follow. See the footnote in Section 4.
  • This is where the present paper and Ref. [8] differ most.

Written By

Paul Papatzacos

Submitted: September 7th, 2020 Reviewed: December 12th, 2020 Published: January 19th, 2021