Porous Flow with Diffuse Interfaces

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.

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 $\nu$ chemical components [2], one obtains

$${\mathcal{F}}^{\mathrm{F}}={\mathcal{F}}^{\mathrm{bF}}+\sum _{\alpha =1}^{\nu }\frac{1}{2}{\mathrm{\Lambda }}^{\alpha }{\left|\mathbf{\text{∇}}{R}^{\alpha }\right|}^2,$$E1

where $R^{\alpha }$ is the upscaled density of fluid component $\alpha$, $\mathcal{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 ${\mathrm{\Lambda }}^{\alpha }$ are constants (Assumption A7, Appendix 1). For later reference ${\mathcal{F}}^{\mathrm{bF}}$ and its differential are

$${\mathcal{F}}^{\mathrm{bF}}=-P^{\mathrm{bF}}+\sum _{\alpha =1}^{\nu }{\mathcal{C}}^{\mathrm{bF}\alpha }{R}^{\alpha },\mathrm{d}{\mathcal{F}}^{\mathrm{bF}}=-{\mathcal{S}}^{\mathrm{bF}}\mathrm{d}T+\sum _{\alpha =1}^{\nu }{\mathcal{C}}^{\mathrm{bF}\alpha }\mathrm{d}{R}^{\alpha },$$E2

where $P^{\mathrm{bF}}$ is the pressure in the bulk fluid, ${\mathcal{S}}^{\mathrm{bF}}$ the entropy per unit volume, and ${\mathcal{C}}^{\mathrm{bF}\alpha }$ the chemical potential of component $\alpha$, divided by its molar mass (the $C^{\alpha }/M^{\alpha }$ 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.¹ 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}^{\circ }$: see Appendix 1, Assumption A0. Refs. [2, 3] assume this limitation.

A note on notation: Right-handed Cartesian coordinates $\left(x_1,x_2,x_3\right)$ are assumed, the plane $\left(x_1,x_2\right)$ being horizontal, and axis $x_3$ pointing upwards. Any vector $\mathbf{A}$ has components $A_k$ where $k$ is 1, 2, or 3. In addition, the notation ${\partial }_t=\partial /\partial t$ and ${\partial }_k=\partial /\partial x_k$ is used. The summation convention applies to latin indexes $i$, $j$, $k$, $l$: an index that is repeated in a term (as in $A_j{B}_j$, or $C_{\mathit{kk}}$) indicates summation (${\sum }_{j=1}^3{A}_j{B}_j$, or ${\sum }_{k=1}^3{C}_{\mathit{kk}}$). Note that an index that is repeated, but not in the same term (as in $A_i=B_i$), means that the expression is valid for all values of the index in the set $\left\{\mathrm{1,2,3}\right\}$. Symbols are otherwise defined when introduced.

2. Pore level equations

The fluid is a mixture of $\nu$ 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:

$${\partial }_t{\rho }^{\alpha }+{\partial }_k\left({\rho }^{\alpha }{v}_k+i_k^{\alpha }\right)=0,\left(\alpha =1,\dots ,\nu \right),$$E3

$${\partial }_t\left(\rho v_i\right)+{\partial }_k\left(\rho v_i{v}_k-t_{\mathit{ki}}\right)=\rho f_i,$$E4

$${\partial }_t\epsilon +{\partial }_k\left(\epsilon v_k+j_k-t_{\mathit{ki}}{v}_i\right)=\rho f_k{v}_k.$$E5

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 ${\rho }^{\alpha }$ is the mass density of component $\alpha$, while $\rho$ is the total density:

$$\rho =\sum _{\alpha =1}^{\nu }{\rho }^{\alpha }.$$E6

Further, $\mathbf{v}$ is the local velocity, defined as the total momentum divided by the total mass, both per volume; ${\mathit{\iota }}^{\alpha }$ is the non-convective mass current of component $\alpha$, with the property

$$\sum _{\alpha =1}^{\nu }{\mathit{\iota }}^{\alpha }=\mathbf{0},$$E7

$\mathit{J}$ is the non-convective energy current, and $t_{\mathit{ij}}$ is the stress tensor; $\epsilon =u+\frac{1}{2}\rho v^2$ is the total energy, $u$ is the internal energy, both per unit volume, and $\mathbf{f}$ is an external force per unit mass of fluid, thus being species-independent. In the case of gravity

$$\mathbf{f}=\mathbf{g},\mathbf{g}=\left(0,0,-g\right),g=9.81\mathrm{m}/{\mathrm{s}}^2.$$E8

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

$$W=-\mathbf{g}\cdot \mathbf{x},\text{whereby}{f}_i=-{\partial }_{i}W.$$E9

The symmetry $t_{\mathit{ij}}=t_{\mathit{ji}}$ 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 $\mathbf{f}$, and the one due the surface stress $\mathbf{n}\cdot \mathbf{t}$ acting on every surface element where the normal vector is $\mathbf{n}$.

Within the van der Waals theory, one expects that ${\iota }_k^{\alpha }$, $J_k$, and $t_{\mathit{ki}}$, 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,

$$t_{\mathit{ij}}=-p{\delta }_{\mathit{ij}}+{\theta }_{\mathit{ij}},$$E10

where $p$ is the pressure and ${\theta }_{\mathit{ij}}$ is the viscous stress tensor (symmetric and linear in the gradients of the $v_i$, 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 $\rho$ 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 ${\partial }_t\rho$ and ${\partial }_k\rho$.

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 $\mathbf{0}=\mathbf{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 $\mathrm{\Sigma }$) and in the solid; likewise, any quantity with superscript S is equal to $0$ on $\mathrm{\Sigma }$ and in the fluid. The quantities in the first line of the table are not defined on $\mathrm{\Sigma }$, but limit values are, as in

Generic	${\rho }^{\alpha }$	$\rho$	$v_k$	${\iota }_k^{\alpha }$	$t_{\mathit{ki}}$	$f_k$	$u$	$j_k$
In fluid	${\rho }^{\mathrm{F}\alpha }$	${\rho }^{\mathrm{F}}$	$v_k^{\mathrm{F}}$	${\iota }_k^{\mathrm{F}\alpha }$	$t_{\mathit{ki}}$	$f_k$	$u^{\mathrm{F}}$	$j_k^{\mathrm{F}}$
In solid		${\rho }^{\mathrm{S}}$	$0$	$0$		$0$	$u^{\mathrm{S}}$	$j_k^{\mathrm{S}}$

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

$${\left({\mathbf{v}}^{\mathrm{F}}\right)}^{\mathrm{\Sigma }}=\mathbf{0},$$E11

expressing that the fluid velocity vanishes on $\mathrm{\Sigma }$. Expressions of the type ${\left(X\right)}^{\mathrm{\Sigma }}$ for some $X$, used in this section and in the next, define the limit of the quantity $X$ at $\mathrm{\Sigma }$, along the line carrying the normal to $\mathrm{\Sigma }$ pointing towards the fluid. The existence of the normal implies some idealisation of the pore surface.

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

$${\rho }^{\mathrm{F}}=\sum _{\alpha =1}^{\nu }{\rho }^{\mathrm{F}\alpha }.$$E12

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$\alpha$ or F, as appropriate, while $B$ is replaced by S. As an example, ${\rho }^{\alpha }\left(\mathbf{x},t\right)$, is a distribution depending on space $\mathbf{x}$ and time $t$. It is continuous in time but not in space, being equal to ${\rho }^{\mathrm{F}\alpha }\left(\mathbf{x},t\right)$ if $\mathbf{x}$ is inside a pore and to ${\rho }^{\mathrm{S}}\left(\mathbf{x},t\right)$ if $\mathbf{x}$ is in the rock. It is not defined when $\mathbf{x}$ is on $\mathrm{\Sigma }$ but we assume that ${\left({\rho }^{\mathrm{F}\alpha }\right)}^{\mathrm{\Sigma }}$ and ${\left({\rho }^{\mathrm{S}}\right)}^{\mathrm{\Sigma }}$ exist.

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

$$\left\{{\partial }_t{\rho }^{\alpha }\right\}+\left\{{\partial }_k\left({\rho }^{\alpha }{v}_k+i_k^{\alpha }\right)\right\}+\left[{\left({\iota }^{\alpha }\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}\right]{\delta }_{\mathrm{\Sigma }}=0,\left(\alpha =1,\dots ,\nu \right).$$E13

From this equation one now gets, using Table 1, the equations that are separately valid inside the pores, on $\mathrm{\Sigma }$, and in the solid:

$$3{\partial }_t{\rho }^{\mathrm{F}\alpha }+{\partial }_k\left({\rho }^{\mathrm{F}\alpha }{v}_k^{\mathrm{F}}+i_k^{\mathrm{F}\alpha }\right)=0,\left(\alpha =1,\dots ,\nu \right),\left(\text{in}F\right)$$E14

$$\left[{\left({\mathit{\iota }}^{\mathrm{F}\alpha }\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}\right]{\delta }_{\mathrm{\Sigma }}=0, \left(\alpha =1,\dots ,\nu \right),\left(\mathrm{on}\mathrm{\Sigma }\right)$$E15

$${\partial }_t{\rho }^{\mathrm{S}}=0.\left(\text{in}S\right)$$E16

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 $\mathbf{0}=\mathbf{0}$. Inside the pores, one simply re-writes Eq. (4), with the superscript F on the density and velocity:

$${\partial }_t\left({\rho }^{\mathrm{F}}{v}_i^{\mathrm{F}}\right)+{\partial }_k\left({\rho }^{\mathrm{F}}{v}_k^{\mathrm{F}}{v}_i^{\mathrm{F}}-t_{\mathit{ki}}\right)={\rho }^{\mathrm{F}}{f}_i,\left(\text{in}F\right).$$E17

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

$$\begin{array}{l}\left\{{\partial }_t\left(u+\frac{1}{2}\rho {\left|\mathbf{v}\right|}^2\right)\right\}+\left\{{\partial }_k\left({\mathit{uv}}_k+\frac{1}{2}\rho {\left|\mathbf{v}\right|}^2{v}_k+j_k-t_{\mathit{ki}}{v}_i\right)\right\}\\ +\left[{\left({\mathit{J}}^{\mathrm{F}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}-{\left({\mathit{J}}^{\mathrm{S}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}\right]{\delta }_{\mathrm{\Sigma }}=f_k\rho v_k.\end{array}$$E18

Using Table 1, one then obtains:

$${\partial }_t\left(u^{\mathrm{F}}+\frac{1}{2}{\rho }^{\mathrm{F}}{\left|{\mathbf{v}}^{\mathrm{F}}\right|}^2\right)+{\partial }_k\left(u^{\mathrm{F}}{v}_k^{\mathrm{F}}+\frac{1}{2}{\rho }^{\mathrm{F}}{\left|{\mathbf{v}}^{\mathrm{F}}\right|}^2{v}_k^{\mathrm{F}}+J_k^{\mathrm{F}}-t_{\mathit{ki}}{v}_i^{\mathrm{F}}\right)=f_k{\rho }^{\mathrm{F}}{v}_k^{\mathrm{F}},\left(\text{in}\mathrm{F}\right)$$E19

$$\left[{\left({\mathit{J}}^{\mathrm{F}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}-{\left({\mathit{J}}^{\mathrm{S}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}\right]{\delta }_{\mathrm{\Sigma }}=0,\left(\mathrm{on}\mathrm{\Sigma }\right)$$E20

$${\partial }_t{u}^{\mathrm{S}}+{\partial }_k{J}_k^{\mathrm{S}}=0,\left(\text{in}\mathrm{S}\right).$$E21

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^{\infty }$ function $m\left(\mathbf{x}\right)$ is introduced, somewhat flat around $\mid \mathbf{x}\mid =0$ and equal to zero for $\mid \mathbf{x}\mid \ge r$, 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, $f\left(\mathbf{x},t\right)$, its average $F\left(\mathbf{x},t\right)$ is obtained by the convolution

$$F\left(\mathbf{x},t\right)\equiv \left(f\ast m\right)\left(\mathbf{x},t\right)={\int }_{{\mathbb{R}}^3}f\left(\mathbf{y},t\right)m\left(\mathbf{x}-\mathbf{y}\right)\mathrm{d}\mathbf{y},$$E22

where the integration is over all of space. The convolution ensures that $F$ is $C^{\infty }$ [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:

$$\left\{{\partial }_t{f}^{\mathrm{Z}}\right\}\ast m={\partial }_t\left(f^{\mathrm{Z}}\ast m\right),$$E23

$$\left\{{\partial }_i{f}^{\mathrm{Z}}\right\}\ast m={\partial }_i\left(f^{\mathrm{Z}}\ast m\right)-{\epsilon }^{\mathrm{Z}}\left[{\left(f^{\mathrm{Z}}\right)}^{\mathrm{\Sigma }}{n}_i{\delta }_{\mathrm{\Sigma }}\right]\ast m,$$E24

where $f^{\mathrm{Z}}\equiv f^{\mathrm{Z}}\left(\mathbf{x},t\right)$ and $\mathrm{Z}$ is F, F$\alpha$, or S; ${\epsilon }^{\mathrm{Z}}$ is 1 if $\mathrm{Z}$ is F or F$\alpha$, −1 if $\mathrm{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:

$${\partial }_t\left({\rho }^{\mathrm{F}\alpha }\ast m\right)+{\partial }_k\left(\left({\rho }^{\mathrm{F}\alpha }{v}_k^{\mathrm{F}}+{\iota }_k^{\mathrm{F}\alpha }\right)\ast m\right)=\left({\left({\mathit{\iota }}^{\mathrm{F}\alpha }\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m,$$E25

$$\left({\left({\mathit{\iota }}^{\mathrm{F}\alpha }\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m=0,$$E26

$${\partial }_t\left({\rho }^{\mathrm{S}}\ast m\right)=0,$$E27

$${\partial }_t\left(\left({\rho }^{\mathrm{F}}{v}_i^{\mathrm{F}}\right)\ast m\right)+{\partial }_k\left(\left({\rho }^{\mathrm{F}}{v}_k^{\mathrm{F}}{v}_i^{\mathrm{F}}-t_{\mathit{ki}}\right)\ast m\right)=-\left({\left(n_k{t}_{\mathit{ki}}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m+f_i\left({\rho }^{\mathrm{F}}\ast m\right).$$E28

Note that the last term on the right-hand side of the last equation originally is $\left({\rho }^{\mathrm{F}}{f}_i\right)\ast m$, but $f_i$ 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 $\mathrm{\Phi }$ is defined first, then follow definitions suggested by the convolution operations in the equations above.

Porosity Let a function $\chi \left(\mathbf{x}\right)$ be 1 when $\mathbf{x}$ is in F, and let it be 0 otherwise.

Then

$$\chi \ast m=\mathrm{\Phi }.$$E29

Note that, according to definition (22), $\mathrm{\Phi }$ can depend on $\mathbf{x}$.

Species adsorption The left-hand side of Eq. (26) defines the amount, $K^{\mathrm{\Sigma }\alpha }$, of component $\alpha$ adsorbed at the pore surface:

$$\left({\left({\mathit{\iota }}^{\mathrm{F}\alpha }\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m=K^{\mathrm{\Sigma }\alpha }.$$E30

Density of solidEq. (27) suggests defining the solid density $R^{\mathrm{S}}$ by:

$${\rho }^{\mathrm{S}}\ast m=\left(1-\mathrm{\Phi }\right)R^{\mathrm{S}}.$$E31

Density of component$\alpha$ The first term on the left-hand side of Eq. (25) suggests defining the density $R^{\alpha }$ of fluid component $\alpha$ by

$${\rho }^{\mathrm{F}\alpha }\ast m=R^{\alpha }.$$E32

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

$$R=\sum _{\alpha =1}^{\nu }{R}^{\alpha }={\rho }^{\mathrm{F}}\ast m.$$E33

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 $V_i$ (without the F superscript since there is no velocity in S or on $\mathrm{\Sigma }$) is then defined by

$$\left({\rho }^{\mathrm{F}}{v}_i^{\mathrm{F}}\right)\ast m=\mathrm{\Phi }{\mathit{RV}}_i.$$E34

Diffusive mass current of component$\alpha$ The second convolution on the left-hand side of Eq. (25) is used to define the upscaled diffusive current in the fluid, denoted $I_k^{\alpha }$ (without the F superscript since there are no diffusive currents in S or on $\mathrm{\Sigma }$). It is not defined as the average of ${\iota }_k^{\mathrm{F}\alpha }$ 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, $I_k^{\alpha }$ is defined by

$$\left({\rho }^{\mathrm{F}\alpha }{v}_k^{\mathrm{F}}+{\iota }_k^{\mathrm{F}\alpha }\right)\ast m=\mathrm{\Phi }{R}^{\alpha }{V}_k+\mathrm{\Phi }{I}_k^{\alpha }.$$E35

Note that, summing this equation over $\alpha$ from 1 to $\nu$, and using previous results, one gets

$$\sum _{\alpha =1}^{\nu }{I}_k^{\alpha }=0.$$E36

Stress tensor The second convolution on the left-hand side of Eq. (28) suggests defining $T_{\mathit{ki}}$, the upscaled version of $t_{\mathit{ki}}$, by

$$\left({\rho }^{\mathrm{F}}{v}_k^{\mathrm{F}}{v}_i^{\mathrm{F}}-t_{\mathit{ki}}\right)\ast m=\mathrm{\Phi }{\mathit{RV}}_i{V}_k-T_{\mathit{ki}}.$$E37

Note that $t_{\mathit{ki}}=t_{\mathit{ik}}$ implies $T_{\mathit{ki}}=T_{\mathit{ik}}$.

Frictional force per unit volume The first convolution on the right-hand side of Eq. (28) suggests defining the frictional force per unit volume ${\mathfrak{F}}_i^{\mathrm{F}}$ by

$$\left({\left(n_k{t}_{\mathit{ki}}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m={\mathfrak{F}}_i^{\mathrm{F}}.$$E38

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:

$${\partial }_t\left[\left(1-\mathrm{\Phi }\right)R^{\mathrm{S}}\right]=0,$$E39

$$K^{\mathrm{\Sigma }\alpha }=0,\left(\alpha =1\dots \nu \right),$$E40

$${\partial }_t\left(\mathrm{\Phi }{R}^{\alpha }\right)+{\partial }_k\left(\mathrm{\Phi }{R}^{\alpha }{V}_k+\mathrm{\Phi }{I}_k^{\alpha }\right)=0,\left(\alpha =1\dots \nu \right),$$E41

$${\partial }_t\left(\mathrm{\Phi }{\mathit{RV}}_i\right)+{\partial }_k\left(\mathrm{\Phi }{\mathit{RV}}_k{V}_i-T_{\mathit{ki}}\right)={\mathfrak{F}}_i^{\mathrm{F}}+\mathrm{\Phi }{\mathit{Rf}}_i.$$E42

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 $\nu$ 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:

$$\begin{array}{l}{\partial }_t\left[\left(u^{\mathrm{F}}+\frac{1}{2}{\rho }^{\mathrm{F}}{\left|{\mathbf{v}}^{\mathrm{F}}\right|}^2\right)\ast m\right]+{\partial }_k\left[\left(u^{\mathrm{F}}{v}_k^{\mathrm{F}}+\frac{1}{2}{\rho }^{\mathrm{F}}{\left|{\mathbf{v}}^{\mathrm{F}}\right|}^2{v}_k^{\mathrm{F}}+J_k^{\mathrm{F}}-t_{\mathit{ki}}{v}_i^{\mathrm{F}}\right)\ast m\right]\\ =\left[{\left({\mathit{J}}^{\mathrm{F}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right]\ast m+f_k\left({\rho }^{\mathrm{F}}{v}_k^{\mathrm{F}}\ast m\right),\end{array}$$E43

$$\left({\left({\mathit{J}}^{\mathrm{F}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m-\left({\left({\mathit{J}}^{\mathrm{S}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m=0,$$E44

$${\partial }_t\left(u^{\mathrm{S}}\ast m\right)+{\partial }_k\left(J_k^{\mathrm{S}}\ast m\right)=-\left[{\left({\mathit{J}}^{\mathrm{S}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right]\ast m.$$E45

As in the case of Eq. (28), Assumption A4 (Appendix 2) has been used to take $f_k$ 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

$$u^{\mathrm{S}}\ast m=\left(1-\mathrm{\Phi }\right)\underset{¯}{{\mathcal{U}}^{\mathrm{S}}},$$E46

Internal energy current of solid

$$j_k^{\mathrm{S}}\ast m=\underset{¯}{J_k^{\mathrm{S}}},$$E47

Solid to fluid energy transfer

$$\left({\left({\mathit{j}}^{\mathrm{S}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m=\underset{¯}{Q_{\mathrm{S}\to \mathrm{F}}},$$E48

Fluid to solid energy transfer

$$\left({\left({\mathit{j}}^{\mathrm{F}}\cdot \mathbf{n}\right)}^{\mathrm{\Sigma }}{\delta }_{\mathrm{\Sigma }}\right)\ast m=-\underset{¯}{Q_{\mathrm{F}\to \mathrm{S}}},$$E49

Internal energy per unit volume of fluid

$$\left(u^{\mathrm{F}}+\frac{1}{2}{\rho }^{\mathrm{F}}{\left|{\mathbf{v}}^{\mathrm{F}}\right|}^2\right)\ast m=\mathrm{\Phi }\underset{¯}{{\mathcal{U}}^{\mathrm{F}}}+\frac{1}{2}\mathrm{\Phi }R{\left|\mathbf{V}\right|}^2,$$E50

Internal energy current of fluid

$$\left(u^{\mathrm{F}}{v}_k^{\mathrm{F}}+\frac{1}{2}{\rho }^{\mathrm{F}}{\left|{\mathbf{v}}^{\mathrm{F}}\right|}^2{v}_k^{\mathrm{F}}+j_k^{\mathrm{F}}-t_{\mathit{ki}}{v}_i^{\mathrm{F}}\right)\ast m=\mathrm{\Phi }{\mathcal{U}}^{\mathrm{F}}{V}_k+\frac{1}{2}\mathrm{\Phi }R{\left|\mathbf{V}\right|}^2{V}_k+\underset{¯}{J_k^{\mathrm{F}}}-T_{\mathit{ki}}{V}_i.$$E51

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

$$\begin{array}{l}{\partial }_t\left(\mathrm{\Phi }{\mathcal{U}}^{\mathrm{F}}+\frac{1}{2}\mathrm{\Phi }R{\left|\mathbf{V}\right|}^2\right)+{\partial }_k\left(\mathrm{\Phi }{\mathcal{U}}^{\mathrm{F}}{V}_k+\frac{1}{2}\mathrm{\Phi }R{\left|\mathbf{V}\right|}^2{V}_k+J_k^{\mathrm{F}}-T_{\mathit{ki}}{V}_i\right)\\ =-Q_{\mathrm{F}\to \mathrm{S}}+\mathrm{\Phi }{\mathit{RV}}_k{f}_k.\left(\text{in}\mathrm{F}\right)\end{array}$$E52

$$Q_{\mathrm{F}\to \mathrm{S}}+Q_{\mathrm{S}\to \mathrm{F}}=0.\left(\mathrm{on}\mathrm{\Sigma }\right)$$E53

$${\partial }_t\left[\left(1-\mathrm{\Phi }\right){\mathcal{U}}^{\mathrm{S}}\right]+{\partial }_k{J}_k^{\mathrm{S}}=-Q_{\mathrm{S}\to \mathrm{F}}.\left(\text{in}\mathrm{S}\right).$$E54

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 $V_i$, summing over $i$, and using Eq. (41). Subtracting the equation thus obtained from Eq. (52), one gets

$${\partial }_t\left(\mathrm{\Phi }{\mathcal{U}}^{\mathrm{F}}\right)+{\partial }_k\left(\mathrm{\Phi }{\mathcal{U}}^{\mathrm{F}}{V}_k+J_k^{\mathrm{F}}\right)=-Q_{\mathrm{F}\to \mathrm{S}}-{\mathfrak{F}}_k^{\mathrm{F}}{V}_k+T_{\mathit{ji}}{\partial }_j{V}_i.$$E55

4. Basis for constitutive equations

In practical application, the upscaled balance equations will be used to calculate, primarily, the densities $R^{\alpha }$, the velocity components $V_i$ 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 $I_k^{\alpha }$, the $T_{\mathit{ij}}$, and the ${\mathfrak{F}}_i$ in terms of $T$ and the $R^{\alpha }$. 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 $T\left(\mathbf{x},t\right)$. Such an equation usually describes the evolution of either total energy or total entropy. In either case, expressions are needed for the currents $J_k^{\mathrm{S}}$ and $J_k^{\mathrm{F}}$. (Expressions for the energy transfers $Q_{\mathrm{S}\to \mathrm{F}}$ and $Q_{\mathrm{F}\to \mathrm{S}}$ 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 $P_{\mathit{ij}}$ 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.²

4.1 The pressure tensor

One considers the upscaled fluid, consisting of a mixture of $\nu$ components in a container with surface $\mathrm{\partial \Omega }$ and volume $\mathrm{\Omega }$, 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 $\mathrm{\partial \Omega }$. The total energy stored in the fluid is

$$F={\int }_{\mathrm{\Omega }}\left[{\mathcal{F}}^{\mathrm{F}}+W\sum _{\alpha =1}^{\nu }{R}^{\alpha }\right]\mathrm{d\Omega },$$E56

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

One now looks for the conditions the $R^{\alpha }$ satisfy when $F$ is at its minimum, given that the total mass, ${\int }_{\mathrm{\Omega }}R\mathrm{d\Omega }$, is constant. This is minimisation with constraint, a standard problem in variational calculus. It is easily found that

$$\frac{\partial {\mathcal{F}}^{\mathrm{bF}}}{\partial R^{\alpha }}+W-\kappa -{\mathrm{\Lambda }}^{\alpha }{\nabla }^2{R}^{\alpha }=0,\left(\text{in}\mathrm{\Omega }\right),$$E57

$$\mathbf{n}\cdot \mathbf{\text{∇}}{R}^{\alpha }=0,\left(\mathrm{on}\mathrm{\partial \Omega }\right),$$E58

where $\kappa$ is a Lagrange multiplier, possibly a function of $T$ but not of the $R^{\alpha }$. 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 ${\partial }_i{R}^{\alpha }$ and summing over $\alpha$. 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:

$${\partial }_k{P}_{\mathit{ik}}-f_i=0.$$E59

where

$$P_{\mathit{ik}}=\left[P^{\mathrm{bF}}-\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }{R}^{\alpha }{\nabla }^2{R}^{\alpha }-\sum _{\alpha =1}^{\nu }\frac{1}{2}{\mathrm{\Lambda }}^{\alpha }{\left|\mathbf{\text{∇}}{R}^{\alpha }\right|}^2\right]{\delta }_{\mathit{ik}}+\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_i{R}^{\alpha }\right)\left({\partial }_k{R}^{\alpha }\right).$$E60

The first equation obviously generalises the classical $\mathbf{\text{∇}}P=f$, and it is natural to generalise Eq. (10) by setting

$$T_{\mathit{ij}}=-\mathrm{\Phi }{P}_{\mathit{ij}}+{\mathrm{\Theta }}_{\mathit{ij}}.$$E61

$T_{\mathit{ij}}$ and $P_{\mathit{ij}}$ are symmetric, so that ${\mathrm{\Theta }}_{\mathit{ij}}={\mathrm{\Theta }}_{\mathit{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 $\mathrm{\Lambda }$‘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 ${\mathrm{\Lambda }}^{\alpha }$ are constants,

$$\mathrm{d}{\mathcal{F}}^{\mathrm{F}}=-{\mathcal{S}}^{\mathrm{bF}}\mathrm{d}T+\sum _{\alpha =1}^{\nu }{\mathcal{C}}^{\mathrm{bF}\alpha }\mathrm{d}{R}^{\alpha }+\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_i{R}^{\alpha }\right)\mathrm{d}\left({\partial }_i{R}^{\alpha }\right).$$E62

This equation leads immediately to two conclusions. There is no additional entropy, and no additional chemical potentials due to large density gradients, since: (i) ${\mathcal{S}}^{\mathrm{F}}$, defined as $\partial {\mathcal{F}}^{\mathrm{F}}/\partial T$, is equal to ${\mathcal{S}}^{\mathrm{bF}}$; (ii) ${\mathcal{C}}^{\mathrm{F}\alpha }$, defined as $\partial {\mathcal{F}}^{\mathrm{F}}/\partial R^{\alpha }$, is equal to ${\mathcal{C}}^{\mathrm{bF}\alpha }$. The differential of ${\mathcal{F}}^{\mathrm{F}}$ can now be re-written with a simpler notation:

$$\mathrm{d}{\mathcal{F}}^{\mathrm{F}}=-{\mathcal{S}}^{\mathrm{F}}\mathrm{d}T+\sum _{\alpha =1}^{\nu }{\mathcal{C}}^{\mathrm{F}\alpha }\mathrm{d}{R}^{\alpha }+\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_i{R}^{\alpha }\right)\mathrm{d}\left({\partial }_i{R}^{\alpha }\right).$$E63

To obtain $\mathrm{d}{\mathcal{S}}^{\mathrm{F}}$ in terms of $\mathrm{d}{\mathcal{U}}^{\mathrm{F}}$ and $\mathrm{d}{\mathcal{F}}^{\mathrm{F}}$ one differentiates ${\mathcal{U}}^{\mathrm{F}}=T{\mathcal{S}}^{\mathrm{F}}+{\mathcal{F}}^{\mathrm{F}}$, obtaining

$$T\mathrm{d}{\mathcal{S}}^{\mathrm{F}}=\mathrm{d}{\mathcal{U}}^{\mathrm{F}}-\sum _{\alpha =1}^{\nu }{\mathcal{C}}^{\mathrm{F}\alpha }\mathrm{d}{R}^{\alpha }-\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_k{R}^{\alpha }\right)\mathrm{d}\left({\partial }_k{R}^{\alpha }\right).$$E64

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

$${\partial }_t\left[\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}+\left(1-\mathrm{\Phi }\right){\mathcal{S}}^{\mathrm{S}}\right]+{\partial }_k\left(\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}{V}_k+P_k\right)=Q,$$E65

where ${\mathcal{S}}^{\mathrm{S}}$ is the entropy of the rock, $P_k$ is a diffusive current, and $Q$ is a source term.

Eq. (64) directly gives the two equations that follow, by replacing $\mathrm{d}$ with ${\partial }_t$ (first equation), then with $V_i{\partial }_i$ (second equation):

$$T{\partial }_t{\mathcal{S}}^{\mathrm{F}}={\partial }_t{\mathcal{U}}^{\mathrm{F}}-\sum _{\alpha =1}^{\nu }{\mathcal{C}}^{\mathrm{F}\alpha }{\partial }_t{R}^{\alpha }-\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_k{R}^{\alpha }\right){\partial }_t\left({\partial }_k{R}^{\alpha }\right).$$E66

$$T{\partial }_i\left({\mathcal{S}}^{\mathrm{F}}{V}_i\right)=T{\mathcal{S}}^{\mathrm{F}}{\partial }_i{V}_i+V_i{\partial }_i{\mathcal{U}}^{\mathrm{F}}-\sum _{\alpha =1}^{\nu }{\mathcal{C}}^{\mathrm{F}\alpha }{V}_i{\partial }_i{R}^{\alpha }-V_i\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_k{R}^{\alpha }\right)\left({\partial }_i{\partial }_k{R}^{\alpha }\right).$$E67

Since it is $\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}$ that is needed to get Eq. (65), one needs to multiply the two equations above with $\mathrm{\Phi }$ and then commute $\mathrm{\Phi }$ with ${\partial }_t$ and ${\partial }_i$. The same commutations are necessary on the right-hand sides of the above equations, so as to get $\mathrm{\Phi }{\mathcal{U}}^{\mathrm{F}}$ 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 ${\partial }_i\mathrm{\Phi }$-terms in the minimal model to be presented. With this simplification in place, and remembering to replace $T_{\mathit{ij}}$ in Eq. (55) with the expression obtained from Eqs. (60) and (61), one obtains, after elementary but somewhat long calculations:

$$\begin{array}{ll}T\left[{\partial }_t\left(\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}\right)+{\partial }_i\left(\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}{V}_i\right)\right]=& -{\partial }_k{J}_k^{\mathrm{F}}+Q_{\mathrm{S}\to \mathrm{F}}-{\mathfrak{F}}_k^{\mathrm{F}}{V}_k+{\mathrm{\Theta }}_{\mathit{ji}}{\partial }_j{V}_i\\ & +\mathrm{\Phi }\sum _{\alpha =1}^{\nu }{\mathcal{C}}^{\mathrm{F}\alpha }{\partial }_i{I}_i^{\alpha }+\mathrm{\Phi }\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_k{R}^{\alpha }\right){\partial }_k{\partial }_i{I}_i^{\alpha }\\ & +\mathrm{\Phi }\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left[R^{\alpha }\left({\partial }_i{V}_i\right){\partial }_k{\partial }_k{R}^{\alpha }+\left({\partial }_k{R}^{\alpha }\right){\partial }_k\left(R^{\alpha }{\partial }_i{V}_i\right)\right].\end{array}$$E68

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:

$$\begin{array}{ll}T\left[{\partial }_t\left(\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}\right)+{\partial }_i\left(\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}{V}_i\right)\right]= & {\partial }_k\left[-J_k^{\mathrm{F}}+\mathrm{\Phi }\sum _{\alpha =1}^{\nu }{A}^{\alpha }{I}_k^{\alpha }\right.\left.+\mathrm{\Phi }\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_k{R}^{\alpha }\right)\left(R^{\alpha }{\partial }_i{V}_i+{\partial }_i{I}_i^{\alpha }\right)\right]\\ & +Q_{\mathrm{S}\to \mathrm{F}}-{\mathfrak{F}}_k^{\mathrm{F}}{V}_k+{\mathrm{\Theta }}_{\mathit{ji}}{\partial }_j{V}_i-\mathrm{\Phi }\sum _{\alpha =1}^{\nu }{I}_i^{\alpha }{\partial }_i{A}^{\alpha },\end{array}$$E69

where $A^{\alpha }={\mathcal{C}}^{\mathrm{F}\alpha }-{\mathrm{\Lambda }}^{\alpha }{\nabla }^2{R}^{\alpha }$. Note that $A^{\alpha }$ occurs twice on the right-hand side of Eq. (69), once as $A^{\alpha }$, and once as ${\partial }_i{A}^{\alpha }$, both times multiplying $I_k^{\alpha }$ in a sum over $\alpha$. Eq. (36) then implies that one can modify the above expression of $A^{\alpha }$ by the addition of any expression that does not depend on $\alpha$. To conform with the notation of Ref. [2], and especially Ref. [3], one then sets

$$A^{\alpha }={\mathcal{C}}^{\mathrm{F}\alpha }+W-{\mathrm{\Lambda }}^{\alpha }{\nabla }^2{R}^{\alpha }.$$E70

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

$$T\mathrm{d}{\mathcal{S}}^{\mathrm{S}}=\mathrm{d}{\mathcal{U}}^{\mathrm{S}},$$E71

which implies

$$T{\partial }_t\left[\left(1-\mathrm{\Phi }\right){\mathcal{S}}^{\mathrm{S}}\right]={\partial }_t\left[\left(1-\mathrm{\Phi }\right){\mathcal{U}}^{\mathrm{S}}\right],$$E72

and, using (54):

$$T{\partial }_t\left[\left(1-\mathrm{\Phi }\right){\mathcal{S}}^{\mathrm{S}}\right]=-{\partial }_k{J}_k^{\mathrm{S}}-Q_{\mathrm{S}\to \mathrm{F}}.$$E73

Taking the sum of Eqs. (69) and (73) one sees that $Q_{\mathrm{S}\to \mathrm{F}}$ cancels, and that the right-hand side consists of a sum of scalar products and of a divergence, say ${\partial }_k{K}_k$. 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 ${\partial }_k{K}_k/T$. This is easily remedied by writing

$$\frac{1}{T}{\partial }_k{K}_k={\partial }_k\left(\frac{K_k}{T}\right)+\frac{K_k}{T^2}{\partial }_{k}T.$$E74

One finally obtains

$${\partial }_t\left[\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}+\left(1-\mathrm{\Phi }\right){\mathcal{S}}^{\mathrm{S}}\right]+{\partial }_k\left(\mathrm{\Phi }{\mathcal{S}}^{\mathrm{F}}{V}_k+\frac{K_k}{T}\right)=Q,$$E75

where

$$K_k=J_k^{\mathrm{F}}+J_k^{\mathrm{S}}-\mathrm{\Phi }{\sum }_{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }\left({\partial }_k{R}^{\alpha }\right)\left(R^{\alpha }{\partial }_i{V}_i+{\partial }_i{I}_i^{\alpha }\right)-\mathrm{\Phi }\sum _{\alpha =1}^{\nu }{A}^{\alpha }{I}_k^{\alpha }.$$E76

$$\mathit{TQ}=-\left(K_k/T\right){\partial }_{k}T+{\mathrm{\Theta }}_{\mathit{ji}}{\partial }_j{V}_i-{\mathfrak{F}}_k^{\mathrm{F}}{V}_k-\mathrm{\Phi }\sum _{\alpha =1}^{\nu }{I}_k^{\alpha }{\partial }_k{A}^{\alpha }.$$E77

$Q_{\mathrm{F}\to \mathrm{S}}$ and $Q_{\mathrm{S}\to \mathrm{F}}$ do not appear in the entropy equation, and are thus not determinable inside the model. Similarly, the currents $J_k^{\mathrm{F}}$ and $J_k^{\mathrm{S}}$ 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. ${\mathfrak{F}}_k^{\mathrm{F}}$ 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 $\mathbf{X}$ (with components $X_i$), and vector currents are denoted $\mathbf{Y}$ (with components $Y_i$), 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 ${\mathrm{\Theta }}_{\mathit{ji}}$ and ${\partial }_j{V}_i$, on the right-hand side of Eq. (77), one sets ${\mathrm{\Delta }}_{\mathit{ji}}\equiv {\partial }_j{V}_i$ and, referring to Appendix 3, especially to Eq. (115), one writes

$${\mathrm{\Delta }}_{\mathit{ij}}=X^A{\delta }_{\mathit{ij}}+X_{\mathit{ij}}^B+{\epsilon }_{\mathit{ijk}}{\tilde{X}}_k,$$E78

$${\mathrm{\Theta }}_{\mathit{ij}}=Y^A{\delta }_{\mathit{ij}}+Y_{\mathit{ij}}^B,$$E79

where $X^A$ and $X_{\mathit{ij}}^B$ are found by replacing $Z$ by $\mathrm{\Delta }$ in Eqs. (115). Note that $Y_{\mathit{ij}}^B$, being symmetric, has no antisymmetric part.

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

$$\begin{array}{l}{X}_k^C=-{\partial }_{k}T\\ X_k^D=-{\mathfrak{F}}_k^{\mathrm{F}}\\ X_k^{\alpha }=-\mathrm{\Phi }{\partial }_k{A}^{\alpha },\end{array}\begin{array}{l}{Y}_k^C=K_k/T\\ Y_k^D=V_k\\ Y_k^{\alpha }=I_k^{\alpha }.\end{array}$$E80

Using Eq. (116), one easily finds that

$$\mathit{TQ}=3X^A{Y}^A+X_{\mathit{ij}}^B{Y}_{\mathit{ij}}^B+{\mathbf{X}}^C\cdot {\mathbf{Y}}^C+{\mathbf{X}}^D\cdot {\mathbf{Y}}^D+\sum _{\alpha =1}^{\nu }{\mathbf{X}}^{\alpha }\cdot {\mathbf{Y}}^{\alpha }.$$E81

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, ${\tilde{X}}_k$, having dropped out.

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

$$\begin{array}{l}{Y}^A=L^{\mathit{AA}}{X}^A+L_{\mathit{kl}}^{\mathit{AB}}{X}_{\mathit{kl}}^B+L_k^{\mathit{AC}}{X}_k^C+L_k^{\mathit{AD}}{X}_k^D+\sum _{\beta =1}^{\nu }{L}_k^{\mathit{A\beta }}{X}_k^{\beta }\\ Y_{\mathit{ij}}^B=L_{\mathit{ij}}^{\mathit{BA}}{X}^A+L_{\mathit{\text{ijkl}}}^{\mathit{BB}}{X}_{\mathit{kl}}^B+L_{\mathit{ijk}}^{\mathit{BC}}{X}_k^C+L_{\mathit{ijk}}^{\mathit{BD}}{X}_k^D+\sum _{\beta =1}^{\nu }{L}_{\mathit{ijk}}^{\mathit{B\beta }}{X}_k^{\beta }\end{array}$$E82

$$\begin{array}{l}{Y}_i^C=L_i^{\mathit{CA}}{X}^A+L_{\mathit{ikl}}^{\mathit{CB}}{X}_{\mathit{kl}}^B+L_{\mathit{ik}}^{\mathit{CC}}{X}_k^C+L_{\mathit{ik}}^{\mathit{CD}}{X}_k^D+\sum _{\beta =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{C\beta }}{X}_k^{\beta }\\ Y_i^D=L_i^{\mathit{DA}}{X}^A+L_{\mathit{ikl}}^{\mathit{DB}}{X}_{\mathit{kl}}^B+L_{\mathit{ik}}^{\mathit{DC}}{X}_k^C+L_{\mathit{ik}}^{\mathit{DD}}{X}_k^D+\sum _{\beta =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{D\beta }}{X}_k^{\beta }\\ Y_i^{\alpha }=L_i^{\mathit{\alpha A}}{X}^A+L_{\mathit{ikl}}^{\mathit{\alpha B}}{X}_{\mathit{kl}}^B+L_{\mathit{ik}}^{\mathit{\alpha C}}{X}_k^C+L_{\mathit{ik}}^{\mathit{\alpha D}}{X}_k^D+\sum _{\beta =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{\alpha \beta }}{X}_k^{\beta },\left(\alpha =1\dots \nu \right).\end{array}$$

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, $L_{\left(n\right)}^{\mathit{MN}}$, obeys

$$L_{\left(n\right)}^{\mathit{MN}}=L_{\left(n\right)}^{\mathit{NM}},\left(M\ne N=A,B,C,D,1,\dots ,\nu \right),$$E83

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 ${\mathrm{\Theta }}_{\mathit{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:

$$L^{\mathit{AZ}}=L^{\mathit{BZ}}=L^{\mathit{ZA}}=L^{\mathit{ZB}}=0,\text{for}\mathrm{any}\mathrm{Z}.$$E84

Then $Y^A=Y_{\mathit{ij}}^B=0$, so that

$${\mathrm{\Theta }}_{\mathit{ij}}=0,$$E85

and the source term of the entropy equation reduces to

$$\mathit{TQ}={\mathbf{X}}^C\cdot {\mathbf{Y}}^C+{\mathbf{X}}^D\cdot {\mathbf{Y}}^D+\sum _{\alpha =1}^{\nu }{\mathbf{X}}^{\alpha }\cdot {\mathbf{Y}}^{\alpha }.$$E86

5.2 Vector currents and forces

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

$$Y_i^C=L_{\mathit{ik}}^{\mathit{CC}}{X}_k^C+L_{\mathit{ik}}^{\mathit{CD}}{X}_k^D+\sum _{\beta =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{C\beta }}{X}_k^{\beta }$$E87

$$Y_i^D=L_{\mathit{ik}}^{\mathit{CD}}{X}_k^C+L_{\mathit{ik}}^{\mathit{DD}}{X}_k^D+\sum _{\beta =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{D\beta }}{X}_k^{\beta }$$E88

$$Y_i^{\alpha }=L_{\mathit{ik}}^{\mathit{C\alpha }}{X}_k^C+L_{\mathit{ik}}^{\mathit{D\alpha }}{X}_k^D+\sum _{\beta =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{\alpha \beta }}{X}_k^{\beta },\left(\alpha =1\dots \nu \right),$$E89

where Onsager symmetry is accounted for. Remembering that $Y_i^{\alpha }=I_k^{\alpha }$, Eq. (36) leads to three constraints on the $L$-coefficients of Eq. (89):

$$\sum _{\alpha =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{C\alpha }}=\sum _{\alpha =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{D\alpha }}=\sum _{\alpha =1}^{\nu }{L}_{\mathit{ik}}^{\mathit{\alpha \beta }}=0.$$E90

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 $L_{\mathit{ik}}^{\mathit{CC}}$ and $L_{\mathit{ik}}^{\mathit{DD}}$.

However, if one limits oneself to the minimal model where: (i) the diffusive entropy current, $K_i/T$, is only due to the ${\partial }_{k}T$ force, (ii) the ${\mathfrak{F}}_k^{\mathrm{F}}$ frictional force is only due to fluid velocity, $V_i$, (iii) each non-convective mass current, $I_i^{\alpha }$, is only due to the gradients of the $A^{\beta }$, then one requires

$$L_{\mathit{ik}}^{\mathit{CD}}=0,\text{and}{L}_{\mathit{ik}}^{\mathit{C\alpha }}=L_{\mathit{ik}}^{\mathit{D\alpha }}=0\text{for}\mathrm{all}\alpha ,$$E91

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 $L_{\mathit{ik}}^{\mathit{\alpha \beta }}$ are zero except when the superscripts are equal: such a matrix would not obey the third constraint in (90). One could of course set $L_{\mathit{ik}}^{\mathit{\alpha \beta }}=0$ for all $\alpha$ and $\beta$, 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, $L_{\mathit{ik}}^{\mathit{CC}}$, $L_{\mathit{ik}}^{\mathit{DD}}$, and $L_{\mathit{ik}}^{\mathit{\alpha \beta }}$, 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:

$$L_{\mathit{ik}}^{\mathit{CC}}={\mathrm{\ell }}^{\mathit{CC}}{\delta }_{\mathit{ik}},L_{\mathit{ik}}^{\mathit{DD}}={\mathrm{\ell }}^{\mathit{DD}}{\delta }_{\mathit{ik}},L_{\mathit{ik}}^{\mathit{\alpha \beta }}={\mathrm{\ell }}^{\mathit{\alpha \beta }}{\delta }_{\mathit{ik}},$$E92

and one gets

$${\mathbf{Y}}^C={\mathrm{\ell }}^{\mathit{CC}}{\mathbf{X}}^C,{\mathbf{X}}^D=\left(1/{\mathrm{\ell }}^{\mathit{DD}}\right){\mathbf{Y}}^D,{\mathbf{Y}}^{\alpha }=\sum _{\beta =1}^{\nu }{\mathrm{\ell }}^{\mathit{\alpha \beta }}{\mathbf{X}}^{\beta }.$$E93

The source term of the entropy equation is now

$$\mathit{TQ}={\mathrm{\ell }}^{\mathit{CC}}{\left|{\mathbf{X}}^C\right|}^2+\left(1/{\mathrm{\ell }}^{\mathit{DD}}\right){\left|{\mathbf{Y}}^D\right|}^2+\sum _{\alpha =1}^{\nu }\sum _{\beta =1}^{\nu }{\mathrm{\ell }}^{\mathit{\alpha \beta }}{\mathbf{X}}^{\alpha }\cdot {\mathbf{X}}^{\beta }.$$E94

Keeping in mind that the ${\mathrm{\ell }}^{\mathit{\alpha \beta }}$-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

$$\sum _{\alpha =1}^{\nu }\sum _{\beta =1}^{\nu }{\mathrm{\ell }}^{\mathit{\alpha \beta }}{\mathbf{X}}^{\alpha }\cdot {\mathbf{X}}^{\beta }=-\sum _{\alpha ,\beta ,\alpha <\beta }{\mathrm{\ell }}^{\mathit{\alpha \beta }}{\left|{\mathbf{X}}^{\alpha }-{\mathbf{X}}^{\beta }\right|}^2.$$E95

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

$${\mathrm{\ell }}^{\mathit{CC}}>0,{\mathrm{\ell }}^{\mathit{DD}}>0,{\mathrm{\ell }}^{\mathit{\alpha \beta }}<0,\left(\alpha <\beta \right).$$E96

These coefficients are determined below.

Using the notation of expressions (80) in the first of expressions (93) one gets $K_i$ proportional to ${\partial }_{i}T$. Now $K_i/T$ is an entropy transport by diffusion, so that $K_i$ is heat transport by diffusion. To recover Fourier’s law one sets

$${\mathrm{\ell }}^{\mathit{CC}}=k/T,\text{implying}{K}_i=-k{\partial }_{i}T,$$E97

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

The second of expressions (93) gives ${\mathfrak{F}}_i^{\mathrm{F}}=V_i/{\mathrm{\ell }}^{\mathit{DD}}$, and it is shown below that the choice

$${\mathrm{\ell }}^{\mathit{DD}}=\frac{K}{{\mathrm{\Phi }}^2\eta },\text{implies}{V}_i=\frac{K}{\mathrm{\Phi }\eta }\left({\partial }_i{P}^{\mathrm{bF}}-{\mathit{Rf}}_i\right),$$E98

away from interphase regions. In this expression, $\eta$ 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 $\eta$ 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:

$$V_i=-\frac{K}{\mathrm{\Phi }\eta }\left[{\partial }_i{P}^{\mathrm{bF}}-{\mathit{Rf}}_i+R\left({\partial }_t{V}_i+V_k{\partial }_k{V}_i\right)-\sum _{\alpha =1}^{\nu }{\mathrm{\Lambda }}^{\alpha }{R}^{\alpha }{\partial }_i{\nabla }^2{R}^{\alpha }\right].$$E99

The sum over $\alpha$ 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 ${\mathit{Rf}}_i$ is of order ${10}^{-3}$, and the material derivative term is of order ${10}^{-9}$. 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:

$$V_i=-\frac{K}{\mathrm{\Phi }\eta }\left[{\partial }_i{P}^{\mathrm{bF}}-{\mathit{Rf}}_i\right]+\mathrm{\Lambda }\frac{\mathit{KR}}{\mathrm{\Phi }\eta }{\partial }_i{\nabla }^{2}R.$$E100

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 $V_i$, that follows by application of the Gibbs-Duhem equation, obtained by differentiating ${\mathcal{F}}^{\mathrm{bF}}$ and using the expression of $\mathrm{d}{\mathcal{F}}^{\mathrm{bF}}$ (Eqs. (2)):

$$V_i=-\frac{K}{\mathrm{\Phi }\eta }\left[\sum _{\alpha =1}^{\nu }{R}^{\alpha }{\partial }_i{A}^{\beta }+{\mathcal{S}}^{\mathrm{F}}{\partial }_{i}T\right],$$E101

where $A^{\alpha }$ 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 $\alpha$:

$$I_k^{\alpha }=-\mathrm{\Phi }\sum _{\beta =1}^{\nu }{\mathrm{\ell }}^{\mathit{\alpha \beta }}{\partial }_k{A}^{\beta },$$E102

where it is reminded that the ${\mathrm{\ell }}^{\mathit{\alpha \beta }}$-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 ${\partial }_k{\mathcal{C}}^{\mathrm{F}\beta }$ behave as $1/R^{\beta }$ when $R^{\beta }$ goes to zero, which is not acceptable when the differential equations of the model are solved numerically. In the minimal model one then sets

$${\mathrm{\ell }}^{\mathit{\alpha \beta }}=-a{\left[R^{\alpha }{R}^{\beta }\right]}^2,\left(\alpha \ne \beta \right),$$E103

using the same positive number $a$ for all elements. Each diagonal element ${\mathrm{\ell }}^{\mathit{\alpha \alpha }}$ is then the negative sum of the off diagonal elements of line $\alpha$. For an example of such a matrix in the case $\nu =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 $\nu$ 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 $\eta$, of Eqs. (97) and (98), are needed as functions of $T$ and of the $R^{\alpha }$. 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 $\alpha$ is (see Refs. [2, 16] and references given there):

$$f^{\mathrm{id}.\mathrm{gas},\alpha }=R_{G}Tln\frac{{\mathcal{V}}^{\alpha }\left(T\right)}{{\mathit{ev}}^{\alpha }},$$E104

where $v^{\alpha }$ is the molar volume, $R_G$ is the gas constant, $e=exp\left(1\right)$, and ${\mathcal{V}}^{\alpha }\left(T\right)$ is a function of temperature with the dimension of a volume per mole. ${\mathcal{V}}^{\alpha }$ 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 ${\mathcal{V}}^{\alpha }$ can also be obtained from an expression of $f^{\mathrm{id}.\mathrm{gas},\alpha }$ in terms of the heat capacity of the component. Assuming ${\mathit{Pv}}^{\alpha }=R_{G}T$, one easily finds that, with $T$ and $v^{\alpha }$ as independent variables, the molar internal energy $u^{\alpha }$ (the “id.gas” superscript is dropped for simplicity) is independent of $v^{\alpha }$, and that one can write

$$\mathrm{d}{u}^{\alpha }=c_V^{\alpha }\left(T\right)\mathrm{d}T,$$E105

where $c_V^{\alpha }$ is the molar heat capacity of component $\alpha$. The molar entropy, $s^{\alpha }$, obeying $T\mathrm{d}{s}^{\alpha }=\mathrm{d}{u}^{\alpha }+P\mathrm{d}{v}^{\alpha }$, one obtains

$$\mathrm{d}{s}^{\alpha }=c_V^{\alpha }\left(T\right)\mathrm{d}T/T+R_G\mathrm{d}{v}^{\alpha }/v^{\alpha }.$$E106

Integrating the last two displayed equations, one obtains $f^{\mathrm{id}.\mathrm{gas},\alpha }$ as $u^{\alpha }-{\mathit{Ts}}^{\alpha }$:

$$f^{\mathrm{id}.\mathrm{gas},\alpha }=u_0^{\alpha }-{\mathit{Ts}}_0^{\alpha }-R_{G}Tln\frac{v^{\alpha }}{v_0^{\alpha }}+{\int }_{T_0}^T\left(1-\frac{T}{T^{\prime }}\right)c_V^{\alpha }\left(T^{\prime }\right)\mathrm{d}{T}^{\prime },$$E107

where $u_0^{\alpha }$ and $s_0^{\alpha }$ are the internal energy and the entropy at a reference state where volume, temperature, and pressure are $v_0^{\alpha }$, $T_0$, and $P_0=R_G{T}_0/v_0^{\alpha }$. Equating the right-hand sides of Eqs. (104) and (107), one obtains an expression for $Tln{\mathcal{V}}^{\alpha }$ in terms of the experimentally measurable function $c_V^{\alpha }$. As stated in Ref. [16], it is in fact the derivative with respect to $T$ of $Tln{\mathcal{V}}^{\alpha }$ that is needed in the differential equations of the thermal model. Assuming $s_0^{\alpha }=0$, one obtains

$$\frac{\mathrm{d}}{\mathrm{d}T}\left[R_{G}Tln{\mathcal{V}}^{\alpha }\left(T\right)\right]=R_{G}ln\frac{{\mathit{eR}}_G{T}_0}{P_0}-{\int }_{T_0}^T\frac{c_V^{\alpha }\left(T^{\prime }\right)}{T^{\prime }}\mathrm{d}{T}^{\prime }.$$E108

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.

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 $f_i$ is approximately constant in the averaging volume

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

A7: The ${\mathrm{\Lambda }}^{\alpha }$ are constant numbers inside the porous medium

A8: The porosity $\mathrm{\Phi }$ 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 $\mathrm{\Sigma }$. For proofs, see [12, 13].

Consider a function $f\left(\mathbf{x},t\right)$ symbolising a physical quantity, continuous in time but dicontinuous in space across a surface $\mathrm{\Sigma }$ that divides space in two regions, called A-side and B-side; $f=f^B$ on the B-side, $f=f^A$ on the A-side. $\mathrm{\Sigma }$ has a normal vector $\mathbf{n}$ at each point, pointing towards A:

$$f\left(\mathbf{x},t\right)=\left\{\begin{array}{ll}{f}^A\left(\mathbf{x},t\right)& \mathbf{x}\mathrm{on}\text{the}\mathrm{A}‐\text{side}\text{of}\mathrm{\Sigma }\\ f^B\left(\mathbf{x},t\right)& \mathbf{x}\mathrm{on}\text{the}\mathrm{B}‐\text{side}\text{of}\mathrm{\Sigma },\end{array}\right.$$E109

One now defines the following regular distributions:

$$\left\{{\partial }_{k}f\left(\mathbf{x},t\right)\right\}=\left\{\begin{array}{ll}{\partial }_k{f}^A\left(\mathbf{x},t\right)& \mathbf{x}\mathrm{on}\text{the}\mathrm{A}‐\text{side}\text{of}\mathrm{\Sigma }\\ {\partial }_k{f}^B\left(\mathbf{x},t\right)& \mathbf{x}\mathrm{on}\text{the}\mathrm{B}‐\text{side}\text{of}\mathrm{\Sigma },\end{array}\right.$$E110

$$\left\{{\partial }_{t}f\left(\mathbf{x},t\right)\right\}=\left\{\begin{array}{ll}{\partial }_t{f}^A\left(\mathbf{x},t\right)& \mathbf{x}\mathrm{on}\text{the}\mathrm{A}‐\text{side}\text{of}\mathrm{\Sigma }\\ {\partial }_t{f}^B\left(\mathbf{x},t\right)& \mathbf{x}\mathrm{on}\text{the}\mathrm{B}‐\text{side}\text{of}\mathrm{\Sigma }.\end{array}\right.$$E111

Then [12, 13].

$${\partial }_{t}f=\left\{{\partial }_{t}f\right\},$$E112

$${\partial }_{k}f=\left\{{\partial }_{k}f\right\}+\left[{\left(f^A\right)}^{\mathrm{\Sigma }}-{\left(f^B\right)}^{\mathrm{\Sigma }}\right]n_k{\delta }_{\mathrm{\Sigma }},$$E113

where ${\delta }_{\mathrm{\Sigma }}$ is a surface-Dirac distribution, its action $⟨{\delta }_{\mathrm{\Sigma }},\phi ⟩$ on a so-called test function $\phi$ [13] being defined as:

$$⟨{\delta }_{\mathrm{\Sigma }},\phi ⟩={\int }_{\mathrm{\Sigma }}\phi \mathrm{d\Sigma }.$$E114

A.3 Appendix 3: Formulas concerning tensors

Let $Z_{\mathit{ij}}$ be an arbitrary second order tensor. It can be shown (see Ref. [15]) that the following expression has general validity:

$$Z_{\mathit{ij}}=Z{\delta }_{\mathit{ij}}+{\widehat{Z}}_{\mathit{ij}}+{\epsilon }_{\mathit{ijk}}{\tilde{Z}}_k,\text{where}\left\{\begin{array}{l}Z=\frac{1}{3}{Z}_{\mathit{kk}},\\ {\widehat{Z}}_{\mathit{ij}}=\frac{1}{2}\left(Z_{\mathit{ij}}+Z_{\mathit{ji}}\right)-\frac{1}{3}{Z}_{\mathit{kk}}{\delta }_{\mathit{ij}},\\ {\tilde{Z}}_k=\frac{1}{2}{\epsilon }_{\mathit{klm}}{Z}_{\mathit{lm}}.\end{array}\right.$$E115

If $Z_{\mathit{ij}}$ is a tensor, then $Z$ is a scalar, $\widehat{Z}$ is a symmetric traceless tensor, and ${\tilde{Z}}_k$ is a pseudo-vector. A most useful property is as follows: if $U_{\mathit{ij}}$ and $V_{\mathit{ij}}$ are two tensors then

$$U_{\mathit{ij}}{V}_{\mathit{ij}}=3\mathit{UV}+{\widehat{U}}_{\mathit{ij}}{\widehat{V}}_{\mathit{ij}}+2{\tilde{U}}_i{\tilde{V}}_i.$$E116