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).
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
The model presented in this chapter was developed in Refs. [1-3]. Ref.  covers the one-component two-phase case, Ref.  is a generalisation to an arbitrary number of components and three phases (two liquids, one gas), Ref.  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. . A partial practical implementation, valid for incomplete wetting, is suggested in Ref.  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 ) 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 , one obtains
where is the upscaled density of fluid component , 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 and its differential are
The van der Waals paper was followed, in 1901, by a paper by Korteweg  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 . 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.  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 . 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.  that the diffuse interface theory presented in Ref.  is consistent for neutral wetting, i.e., for pore level wetting angles around : see Appendix 1, Assumption A0. Refs. [2, 3] assume this limitation.
A note on notation: Right-handed Cartesian coordinates are assumed, the plane being horizontal, and axis pointing upwards. Any vector has components where is 1, 2, or 3. In addition, the notation and is used. The summation convention applies to latin indexes , , , : an index that is repeated in a term (as in , or ) indicates summation (, or ). Note that an index that is repeated, but not in the same term (as in ), means that the expression is valid for all values of the index in the set . 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 . Greek superscripts indicate the species so that is the mass density of component , while is the total density:
Further, 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
is the non-convective energy current, and is the stress tensor; is the total energy, is the internal energy, both per unit volume, and 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 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 , and the one due the surface stress acting on every surface element where the normal vector is .
Within the van der Waals theory, one expects that , , and , 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 is the pressure and is the viscous stress tensor (symmetric and linear in the gradients of the , 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 . 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 and .
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 -values indicate no material transport in the solid, the third zero value follows from Assumption A4 (Appendix 1): as shown by Marle , this assumption implies that the momentum balance equations in the solid and on the pore surface have the form . It is important to keep in mind, especially when averaging, that any quantity with an F superscript is equal to on the pore surface (denoted by ) and in the solid; likewise, any quantity with superscript S is equal to on and in the fluid. The quantities in the first line of the table are not defined on , but limit values are, as in
expressing that the fluid velocity vanishes on . Expressions of the type for some , used in this section and in the next, define the limit of the quantity 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 is replaced by For F, as appropriate, while is replaced by S. As an example, , is a distribution depending on space and time . It is continuous in time but not in space, being equal to if is inside a pore and to if is in the rock. It is not defined when is on but we assume that and exist.
From this equation one now gets, using Table 1, the equations that are separately valid inside the pores, on , and in the solid:
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 . Inside the pores, one simply re-writes Eq. (4), with the superscript F on the density and velocity:
Using Table 1, one then obtains:
The Marle averaging process  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 , large when compared to a pore radius, small when compared to a linear dimension of the reservoir. A function is introduced, somewhat flat around and equal to zero for , normalised so that its integral over all space is equal to . (See Eq. (17) in Ref.  for an example of such a function.) Given any function of space and time, , its average is obtained by the convolution
where the integration is over all of space. The convolution ensures that is .
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 . Step 2: the following rules  are applied, allowing to take the differential operators out of the averaging convolutions:
where and is F, F, or S; is 1 if is F or F, −1 if 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.
Note that the last term on the right-hand side of the last equation originally is , but 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 be 1 when is in F, and let it be 0 otherwise.
Note that, according to definition (22), can depend on .
Species adsorption The left-hand side of Eq. (26) defines the amount, , of component adsorbed at the pore surface:
Density of solidEq. (27) suggests defining the solid density by:
Density of componentThe first term on the left-hand side of Eq. (25) suggests defining the density of fluid component by
Note that Eq. (12) implies that the averaged total fluid density, , 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 (without the F superscript since there is no velocity in S or on ) is then defined by
Diffusive mass current of componentThe second convolution on the left-hand side of Eq. (25) is used to define the upscaled diffusive current in the fluid, denoted (without the F superscript since there are no diffusive currents in S or on ). It is not defined as the average of because both pore level effects of convection and diffusion contribute to it . Keeping in mind the constraint that the averaged equations should have the same form as the original ones, 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 , the upscaled version of , by
Note that implies .
Frictional force per unit volume The first convolution on the right-hand side of Eq. (28) suggests defining the frictional force per unit volume 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).
As in the case of Eq. (28), Assumption A4 (Appendix 2) has been used to take 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
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 , summing over , 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 , the velocity components and the temperature . Assuming that the approximation of constant 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 , the , and the in terms of and the . 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.  and Appendix 1), from which it follows that just one additional equation is needed for calculating . Such an equation usually describes the evolution of either total energy or total entropy. In either case, expressions are needed for the currents and . (Expressions for the energy transfers and 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 that replaces the usual scalar pressure.
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 . 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
One now looks for the conditions the satisfy when is at its minimum, given that the total mass, , 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 but not of the . See Ref.  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 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 , and it is natural to generalise Eq. (10) by setting
and are symmetric, so that . 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
This equation leads immediately to two conclusions. There is no additional entropy, and no additional chemical potentials due to large density gradients, since: (i) , defined as , is equal to ; (ii) , defined as , is equal to . The differential of can now be re-written with a simpler notation:
To obtain in terms of and one differentiates , 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 is the entropy of the rock, is a diffusive current, and is a source term.
Eq. (64) directly gives the two equations that follow, by replacing with (first equation), then with (second equation):
Since it is that is needed to get Eq. (65), one needs to multiply the two equations above with and then commute with and . The same commutations are necessary on the right-hand sides of the above equations, so as to get 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 -terms in the minimal model to be presented. With this simplification in place, and remembering to replace 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:
where . Note that occurs twice on the right-hand side of Eq. (69), once as , and once as , both times multiplying in a sum over . Eq. (36) then implies that one can modify the above expression of by the addition of any expression that does not depend on . To conform with the notation of Ref. , and especially Ref. , 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
and, using (54):
Taking the sum of Eqs. (69) and (73) one sees that cancels, and that the right-hand side consists of a sum of scalar products and of a divergence, say . A last step remains because the left-hand side has a multiplicative factor (see Eqs. (69) and (73)). Dividing both sides by does not affect the scalar products but it modifies the divergence, at least if one assumes that is variable, to . This is easily remedied by writing
One finally obtains
and do not appear in the entropy equation, and are thus not determinable inside the model. Similarly, the currents and 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. 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 (with components ), and vector currents are denoted (with components ), with a superscript to discriminate between the different currents and forces in an alphabetic order aa shown below. (has been chosen instead of the traditional 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 and , on the right-hand side of Eq. (77), one sets and, referring to Appendix 3, especially to Eq. (115), one writes
where and are found by replacing by in Eqs. (115). Note that , 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 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, , 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 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 subscripts, denoted say, , obeys
where the 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
With hindsight, one knows that the upscaled viscosity tensor is not required, so that one is justified in setting equal to zero all the ‘s in the first two lines of the system of Eqs. (82), and also setting to zero all the ‘s related to the zeroed ones by the Onsager symmetry:
Then , 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:
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 coefficients. Also, thermal conduction and permeability can be modelled by second order tensors, through tensors and .
However, if one limits oneself to the minimal model where: (i) the diffusive entropy current, , is only due to the force, (ii) the frictional force is only due to fluid velocity, , (iii) each non-convective mass current, , is only due to the gradients of the , 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 are zero except when the superscripts are equal: such a matrix would not obey the third constraint in (90). One could of course set 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, , , and , 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. ). 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 proportional to . Now is an entropy transport by diffusion, so that is heat transport by diffusion. To recover Fourier’s law one sets
where is the thermal conductivity of the averaged solid–fluid system. See Ref. .
The second of expressions (93) gives , and it is shown below that the choice
away from interphase regions. In this expression, is the viscosity of the pore fluid, and 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. ). The possibility of letting 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. . 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 is of order , and the material derivative term is of order . 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.  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.  that relative permeabilities cannot capture the full complexity of two-phase flow.
Given below is another version of , that follows by application of the Gibbs-Duhem equation, obtained by differentiating and using the expression of (Eqs. (2)):
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 behave as when goes to zero, which is not acceptable when the differential equations of the model are solved numerically. In the minimal model one then sets
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. .
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.  for a detailed presentation.
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.  and references given there).
It is shown in Ref.  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 that drops out under differentiation when 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 is the molar volume, is the gas constant, , and is a function of temperature with the dimension of a volume per mole. depends on the atomic masses, the principal moments of inertia of the molecule, .... It can be obtained from statistical physics, as shown in Ref. .
It is briefly shown below that can also be obtained from an expression of in terms of the heat capacity of the component. Assuming , one easily finds that, with and as independent variables, the molar internal energy (the “id.gas” superscript is dropped for simplicity) is independent of , and that one can write
where is the molar heat capacity of component . The molar entropy, , obeying , one obtains
Integrating the last two displayed equations, one obtains as :
where and are the internal energy and the entropy at a reference state where volume, temperature, and pressure are , , and . Equating the right-hand sides of Eqs. (104) and (107), one obtains an expression for in terms of the experimentally measurable function . As stated in Ref. , it is in fact the derivative with respect to of that is needed in the differential equations of the thermal model. Assuming , one obtains
Concerning examples of numerical solutions of the equations of the minimal model, see Ref.  for phase segregation, and for coning at uniform temperature; see Ref.  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 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
Consider a function 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; on the B-side, on the A-side. has a normal vector at each point, pointing towards A:
One now defines the following regular distributions:
where is a surface-Dirac distribution, its action on a so-called test function  being defined as:
A.3 Appendix 3: Formulas concerning tensors
Let be an arbitrary second order tensor. It can be shown (see Ref. ) that the following expression has general validity:
If is a tensor, then is a scalar, is a symmetric traceless tensor, and is a pseudo-vector. A most useful property is as follows: if and are two tensors then
- 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.  differ most.