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

## Abstract

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.

### Keywords

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

where

where

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

**A note on notation:** Right-handed Cartesian coordinates

## 2. Pore level equations

The fluid is a mixture of

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

Further,

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

The symmetry

Within the van der Waals theory, one expects that

where

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

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

Generic | ||||||||
---|---|---|---|---|---|---|---|---|

In fluid | ||||||||

In solid |

expressing that the fluid velocity vanishes on

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

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

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

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

where the integration is over all of space. The convolution ensures that

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

where

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

**Porosity** Let a function

Then

Note that, according to definition (22),

**Species adsorption** The left-hand side of Eq. (26) defines the amount,

**Density of solid**Eq. (27) suggests defining the solid density

**Density of component**

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

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

**Diffusive mass current of component**

Note that, summing this equation over

**Stress tensor** The second convolution on the left-hand side of Eq. (28) suggests defining

Note that

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

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

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

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

## 4. Basis for constitutive equations

In practical application, the upscaled balance equations will be used to calculate, primarily, the densities

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 use* expression (1) *in both derivations.* ^{2}

### 4.1 The pressure tensor

One considers the upscaled fluid, consisting of a mixture of

where

One now looks for the conditions the

where

The next step consists in multiplying Eq. (57) with

where

The first equation obviously generalises the classical

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

This equation leads immediately to two conclusions. There is no additional entropy, and no additional chemical potentials due to large density gradients, since: (i)

To obtain

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

Eq. (64) directly gives the two equations that follow, by replacing

Since it is

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

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

One finally obtains

where

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

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

where

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

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

In these expressions, the

where the

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

Then

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

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

However, if one limits oneself to the minimal model where: (i) the diffusive entropy current,

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

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 one gets

The source term of the entropy equation is now

Keeping in mind that the

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

where

The second of expressions (93) gives

away from interphase regions. In this expression,

The sum over

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

where

where it is reminded that the

using the same positive number

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

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,

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

where

It is briefly shown below that

where

Integrating the last two displayed equations, one obtains

where

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

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

**A7:** The

**A8:** The porosity

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

Consider a function

One now defines the following regular distributions:

where

## A.3 Appendix 3: Formulas concerning tensors

Let

If

## Notes

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