Non-isothermal transport of multi-phase fluids in porous media. Constitutive equations

We develop constitutive equations for multi-component, multi-phase, macro-scale flow in a porous medium exposed to temperature-, composition-, and pressure -gradients. The porous medium is non-deformable. We define the pressure and the composition of the representative elementary volume (REV) in terms of the volume and surface averaged pressure and the saturation, and the respective driving forces from these variables. New contributions due to varying porosity or surface tension offer explanations for non-Darcy behavior. The interaction of a thermal and mechanical driving forces give thermal osmosis. An experimental program is suggested to verify Onsager symmetry in the transport coefficients.


Introduction
We have recently [1] derived the entropy production, σ, of a representative elementary volume (REV) in a heterogeneous system, a multi-phase, multi-component fluid in a porous non-deformable matrix. The coarse-grained description of the REV was formulated for systems that obey Euler homogeneity. A Gibbs equation could therefore be formulated for the REV itself, and used as a starting point, as is common in non-equilibrium thermodynamics [2]. Once the entropy production was found from this and the balance equations, the driving forces and the constitutive equations can be found. This will be done here.
We shall see that the form of the equations is the same as for homogeneous systems, but the driving forces are particular for the multi-component, non-isothermal and multi-phase flow in porous media. We are seeking internal area to volume (no film formation). For such systems we proceed to find detailed expressions for the part of the chemical potential of the REV. The driving forces, due to temperature -, pressure -and chemical potential gradients, are specified. The simple case of flow of two immiscible fluids is used to bring out the specifics of porous media.
An experimental program is suggested in the end to verify Onsager symmetry in the transport coefficients.

Thermodynamic variables for the REV
The central concept in this analysis is the representative elementary volume; the REV [12,13]. Its characteristic size, l, is small compared to the size (length) of the full system, L, but large compared to the characteristic pore length and diameter. An illustration of the REV is given by the squares in Fig. 1. The REV (square) consists of several phases and components. The problem is to find the representation on the larger scale. For each point in the porous system, represented by the (blue) dot in Fig.1, we use the REV around it to obtain the variables (U REV , S REV , V REV , M REV i ) of the REV. This was first defined in [1]. From the Euler homogeneity of these variables, the possibility followed to define also the temperature, pressure and chemical potentials of the REV, (T, p, µ i ) [1]. Figure 1: A representative elementary volume (REV)(magenta squares, length l) in a system (white box of length L) that is exposed to for instance a pressure difference, ∆p, a temperature difference, ∆T , and differences in chemical potentials ∆µ i . A (blue) dot is used to represent the state that characterizes the REV. The set of state variables are positioned on the x-axis: The temperature, pressure and the chemical potentials, (T, p, µ i )(x).
The procedure will be recapitulated. We follow standard thermodynamics and choose as a basis set of variables, the set that consists of internal energy, entropy, volume, (U, S, V ), and the component masses, M i . These variables, are given superscript REV, and constitute the only independent variables of the REV.
The value of each of these variables of the REV is obtained as a sum of contributions from each phase, interface and three-phase contact line present [12,13]. The contributions arise say from the phase volumes, the surface areas and the contact line lenghts and are pore-scale variables; they are not independent variables on the macro-scale.
Assuming Euler homogeneity, means that a REV of the double size, for example, has double the energy, entropy, and mass, as well as double the surface areas of various types and double the line lengths. The average surface area, pore length and pore radius, as well as the curvature of the surfaces in the REV, are then everywhere the same.
A system of k components in m homogeneous phases, has a volume, V REV , with contributions from the homogeneous bulk phases V α,REV , m ≥ α ≥ 1, and the excess line volumes, V αβδ,REV , m ≥ α > β > δ ≥ 1.
The excess surface volumes are by construction zero. The excess line volumes are not, because the dividing surfaces in general cross each other along three different lines. We neglect these contributions, which are normally small also in porous media. The volume of the pores is Superscript p is used for pore. The porosity, φ, and the degree of saturation,Ŝ α (saturation for short), are Superscript α is used for a component, which is equal to the phase in the present case. The porosity and the saturation do not depend on the size of the REV, and have therefore no REV superscript.
We shall often use the example of two immiscible one-component phases w and n in a solid porous material r of porosity φ, where contact line contributions are negligible. We can think of phase w as wetting, and n as non-wetting. The mass variables are then from Eq. 4: When an interface is formed between two phases, we are free to choose the position of the interface such that one of the components has a zero excess mass. This is the position of the equimolar surface of this component. This position is convenient because the number of variables are reduced. When we use the equimolar surface of n, M rn n = M wn n = 0, and when we use the equimolar surface of w at the surface of the solid, M rw w = 0. These choices simplify the description of the REV.
The expressions for U REV and S REV are similar to Eq. 4. This way to construct a REV is reminiscent of the geometric construction of a state function, proposed for flow in porous media by McClure et al. [14].
The basis set of macro-scale variables of the REV (U REV , S REV , V REV , M REV i ) apply to the whole REV. The temperature, pressure and chemical potentials of the REV, (T, p, µ i ), on the macro-scale were next defined, as is normal in thermodynamics, as partial derivatives of the internal energy. These definitions are normal in the sense that they have the same form as they have in homogeneous systems. They differ from definitions in normal homogeneous systems in that the variables here (say U REV ) have contributions from all parts of the heterogeneous REV. The intensive variables T, p and µ i are not averages of the corresponding variables on the pore-scale, as was also pointed out by Hassanizadeh and Gray [11,12].
The macro-scale densities of internal energy, entropy and mass; in the example, u, s, ρ i , do not depend on the size of the REV. The densities are therefore convenient when we need to integrate across the system [1]. They are, however, functions of the position of the REV.

The entropy production of two-phase flow
Pressure-driven mass flow through porous media can lead to gradients in composition and temperature, and vice versa; temperature gradients can lead to mass flow, separation and pressure gradients. The interaction of such flows is of interest, and motivated the search for convenient forms of the entropy production [1].

Component flow variables
From the Gibbs equation for the REV, we derived the entropy production for transport of heat and two immiscible fluid phases through the REV [1]. With transport in the x-direction only, the entropy production σ of the example system was The frame of reference for the mass transport is the non-deformable solid matrix, J r ≡ 0. Here J q is the sensible heat flux in J/(m 2 s)), T is the temperature (in K), J i is a component flux (in kg.m −2 .s −1 ) and ∂µ i,T /∂x is the gradient of chemical potential (in J.kg −1 .m −1 ) evaluated at constant temperature. All properties are for the REV so the superscript is omitted.
The thermal force conjugate to the heat flux is the gradient of the inverse temperature, where the temperature was defined for the REV as a whole, see [1]. This force will not be further discussed.
The chemical force conjugate to the mass flux is the inverse temperature times the negative gradient of µ i,T .
This driving force has several contributions, which must be defined for the REV. We are seeking more specific expressions for the last driving forces in Eq.6.
In order to take all REV variables into account, see Section 4 below, we derive the chemical potential starting from the Gibbs energy, G: where G applies to the REV and G i is defined for component i in the last identity. The total differential of U is introduced and we obtain The full chemical potential is the derivative with respect to M i , which was defined in the previous section; The total differential of the chemical potential is: T,M k are partial specific quantities. The last term describes the change in the chemical potential by changing composition of the medium. The dµ i,T is now defined as a part of the whole differential: The last term is zero when the composition is uniform. The pressure gradient can be introduced as a driving force in Eq. 6 through this equation.

Volume flow as variable
The volume flow is often measured, and is thus a central variable. In the simple case of a single fluid, say w, we obtain from Eq.6 and 11 that The volume flow is J V = J w V w for a single fluid. For two fluids In the absence of a gradient in composition, dµ c i = 0, the same expression applies with this volume flow. We may also follow Hansen et al [15] and write the component contributions as where the saturation has been introduced, and v i is the volume flow of i.
With two fluids in a uniform, non-deformable rock, there are three components. On the coarse-grained level, these are mixed. We assume that dµ c r = 0, and obtain Gibbs-Duhem's equation on the form where ρ i is the density of i in the REV (in kg.m 3 ).
This can be used with Eq.11 and J V to change Eq.6 into The entropy production is invariant, and this defines J D as the velocity of component w relative to n (in m.s −1 ): The entropy production 15 has also three terms. While the first term on the right-hand side is the same as before, the second term is the volume flow with minus the pressure gradient over the temperature as driving force, and the third term is the velocity difference with the chemical potential gradient times the density over the temperature as driving force.
Equations 6 and 15 are equivalent. They describe the same entropy production or flow dissipation. They provide alternative choices of conjugate thermodynamic force-flux pairs. The choice to use in the particular case, is determined by practical reasons; what can be measured or not, which terms are zero. For instance, under isothermal conditions we need not take the term containing the heat flux along, even if heat may be transported reversibly. One set may give a negative contribution to the entropy production (work is done), but the overall entropy production is positive, of course. Each set can be used to obtain constitutive equations for transport on the macro-scale. We shall proceed to find these for porous media flow, finding first more detailed expressions for the driving forces. In order to do so, we again use the additive properties of the REV-variables. 4 The chemical driving force

Saturation-dependent contributions
We need the specific contribution to the chemical potential gradients in the entropy production in Eqs. 6 The expression gives the Gibbs energy contributions of component i to the REV. We neglected again possible contributions from contact lines.
In the case of two immiscible, one-component fluids in a non-deformable porous rock, we obtain for the nonwetting fluid.
For immiscible, one-component fluids the label indicating the components also gives the phase. The densities g n n , g rn n and g wn n are averages over V n,REV , Ω rn,REV and Ω wn,REV . Their local values in the pores may vary around these averages.
The concentration dependent part of the chemical potential of i, µ i , is (in J.kg −1 ) Here R is the gas constant (in J.K −1 .mol −1 ) and W i is the molar mass (in kg.mol −1 ). The chemical potential is measured referred to a standard state, µ 0 i , having the local concentration ρ 0 i in all the pores. In the description of porous media, a convenient reference state may be the state when one component is filling all pores, or the saturation is unity,Ŝ i = 1. The mass density of i in the REV for the standard state is ρ 0,REV By introducing these definitions into Eq.19, we obtain for the concentration dependent part of the driving force We see that any variation in saturation between REVs along the system, will lead to a driving force. We integrate between two REVs, and we obtain the chemical driving force for porous media flow We have seen that the return to the Gibbs energy of the porous medium helped define the chemical potential in terms of properties relevant to porous media. All variables are measurable.

The pressure of a REV
We find the contribution from the pressure to the driving force, by starting as above with the extensive property that holds the variable. For the pressure, this is the grand potential. The compressional energy of the REV is equal to minus the grand potential: The grand potential of the REV is additive, which gives We introduce contributions from all phases, surfaces and lines. This gives for the pressure of the REV: The superscript denotes the relevant phase. The last equation makes it possible to compute the pressure of the REV, p, from the pressures in the bulk phases, and the surface-and line tensions. With knowledge of the pressure in the REV, we find the driving force, −dp/dx, in the entropy production, Eq.15. We explain now how we can define and compute the pressure from Eq.25, using the example of two immiscible fluid in a non-deformable medium.
We follow Eq.25 and sum over the n, w and r bulk phases, and the nr, wr and nw−interfaces. The situation can be illustrated for a single cylindrical pore, see The pressure of the REV from Eq. 25 can then be written as: Contact-line contributions were again not taken along, for simplicity. A consequence of the porous medium being homogeneous is that Ω rp,REV /V REV is the same everywhere. The ratio can be used as a measure of the average curvature of the pore surface, as will be explained below.
The volume-averaged contributions to the pressure from the homogeneous phases is given the symbolp: The first term in the last equality shows that the saturation gives an important contribution to the volume-averaged pressure. The contributions of the 2nd and 3rd terms are due to p w and p r . These contributions are usually constant.
The surface-averaged contributions to the pressure are likewise given a separate symbol: The contribution ofp c to the pressure, p, is often called the capillary pressure. With an (approximately) cylindrical pore geometry, we can define the average radius of the pores by This may be a good assumption in the absence of film formations. By introducing r into Eq.29 we obtain the capillary pressurep Again the first term shows that saturation gives an important contribution. The 2nd term only depends on the temperature and is usually constant. The 3rd term is proportional to the surface area of the fluid-fluid interface. In many experiments this surface area is much smaller than Ω rp,REV . When that is the case, this term is negligible.
The effective pressure of the REV is thus, for short: The three equations above give an expression for the REV pressure p for the example system.
Approximations should be tested with the more detailed expressions. To estimate the size of the various contributions, it is convenient to use mechanical equilibrium for the contact line and for the surface, although this condition may not apply to the REV, not even under stationary flow conditions. With a balance of forces at the three-phase contact lines, Young's law applies for the surface tensions: γ nr − γ wr = γ nw cos θ, where θ is the (average) contact angle. When there is furthermore mechanical equilibrium at the fluid-fluid interfaces, the pressure difference between the fluid is given by Young-Laplace's law, p n − p w = 2γ nw cos θ/r.
In the single-fluid (w) case, Eq.27 simplifies. The volume-averaged contribution becomes, and by introducing Eq.31, we obtain In this section, we have defined in detail what we mean by the pressure of a REV. We have found, using the grand potential, that it can be regarded as result of volume-and surface average properties. These contributions enter the driving force in Eqs. 6 and 15, to be further discussed below.

The pressure difference as driving force across a porous medium
The driving force for volume flow is the gradient of the REV pressure. To measure the pressure p inside the REV is difficult. The pressure in the fluid phases adjacent to the porous medium can be determined. Tallakstad et al. [16] defined the measured pressure difference, ∆p , at steady state, as an average over the value ∆p(t) over the time of measurement: Here t is the time and ∆ refers to the extension of the system. Subscript 's' denotes the start and 'e' denotes the end of the measurement. We will take this pressure difference as our ∆p.
The pressure differences ∆p w and ∆p n can also be found when there is continuity in the fluids, w and n, respectively.
By taking the difference between the inlet and outlet in Eq.32, we have an interpretation of the pressure difference; The question is now how we can assess the right-hand side of this equation, using Eqs. 28 and 29.

Large pressure differences
When the pressure drop across the porous plug is large compared to the capillary pressure contribution, the surface contributions and thereforep c can be neglected. Furthermore p n = p w . In the pressure difference, the terms with constant φ and p r disappears, and the pressure difference is: The pressure is applied to the whole cross-sectional area. This explains that the net driving force becomes a fraction, φ, of ∆p w . In other words, the force applies to the fraction φ of the pore area.
The conditions leading to Eq.37 are common in the laboratory. Some numerical values for the air-glycerol system, [4], can illustrate when the conditions apply. The value of 2φγ wr /r is of the same order of magnitude as p cŜn φ (400 Pa) when the surface tension γ = 6.4 · 10 −2 N m −2 , the average pore radiusr = 0.2 mm and the porosity φ = 0.63. A typical value of ∆p in the experiments is close to 30 kPa, which is far from the limit where capillary effects are significant.

Small pressure differences
For small capillary numbers the effective pressure drop across a porous plug is comparable to or smaller than the capillary pressure. Surface contributions need be taken into account. Equation 36 gives the effective pressure difference. When we can assume a constant average radiusr, and constant porosity, we obtain A fluid will be transported when the surface tensions of the fluids with the wall are different and there is a difference in the saturation. When there is only one fluid in the porous medium, cf. Eq.34, and we have constantr and porosity, the pressure difference becomes The last term can lead to mass transport, when the surface tension changes.

Constitutive equations
The constitutive equations follow from the entropy production. We present these on differential form for two incompressible flows. From Eq. 6 we have: We can also use Eq. 15 and obtain The flux-force relations are linear on this level. The two conductivity matrices can be expressed in each other. The element l qq is the same in both formulations. When the REV can be regarded as a thermodynamic state [4,5], we can expect that the conductivity matrices on this local form are symmetric, or that the Onsager relations apply, l ij = l ji . FlekkÃ¸y et al. [19] and Pride et al. [20] argued that the Onsager relations are obeyed also on the REV level. Experimental proof for the Onsager relations of two-phase flow in porous media does not yet exist, however, and we shall see below how this possibly can be achieved.
We have discussed above how the overall driving forces can be determined. We need to integrate across the system, in order to study their effect on experiments. We integrate the linear laws 41 across the REV, and obtain and Here L ij ≡ l ij /l and l is the length of the REV, and the driving forces are defined by Eqs. 36 and 22.
The coefficients may become dependent on the force through the integration. To illustrate this, consider an example. Two fluids in a capillary were studied, and linear laws were first written for each of them on the pore-level [21].
The average velocity of a bubble was found by integration from the pore-to the macro-level, using the configurational distribution Π(x b ), of the position, x b , of the fluid interfaces. The capillary pressure depended on x b . The averaging procedure gave the conductivity as a function of (∆p -∆p c ) in the terminology of this paper.
In the remainder of the paper we will discuss experimental conditions that allow us to determine these coefficients.
The presentation follows closely the derivation of Stavermann [22] and Katchalsky and coworkers [23] for transport in discrete systems of polymer membranes, see also F/orland [24]. We refer to these works for further definitions of transport coefficients.

Isothermal, single fluid flow.
For an isothermal single fluid w, flowing inside a porous medium, the entropy production 43 has one term; the volume flow times the negative pressure difference over the temperature. By including the constant temperature in the transport coefficient, we obtain the common linear law. With the permeability L p , we write where L p ≡ L pp /T . The permeability is normally a function of state variables (pressure, temperature). In the hydrodynamic regime it is a function of viscosity, L p = L p (p, T, η). By introducing the new expression for the pressure 34, we obtain When the permeability and porosity are constant, L V V ≡ L p φ. The equation predicts a threshold value for flow if there is a (significant) change in the surface tension across the REV. Transport will take place, when ∆p w > 2∆γ wr /r. The permeability L V V is inversely proportional to the viscosity η of the fluid in the hydrodynamic regime. Interestingly, Boersma et al. [8] and Miller et al. [7] plotted the volume flow versus the hydrostatic pressure difference ∆p w and found a deviation from Darcy's law in the form of a pressure threshold, for water or water solutions in clay. They offered no explanation for this. Also Bernadiner et al. [9] and Swartzendruber [6] plotted the volume flow of water solution J V versus the pressure gradient in sandstone with low clay content [9], and in NaClsaturated Utah bentonite [6]. The thresholds that they observed depended on the content of salt in the permeating solution. They explained the thresholds by water adsorption and pore clogging by colloids [9]. According to Eq.45, a varying surface tension (due to a varying adsorption and clogging) might explain the existence of a threshold or a non-linear flux-force relation. The dependence of the coefficient L V V on the threshold pressure can also have other explanations, cf. the case described above [21]. This non-linearity does not prevent the use of non-equilibrium thermodynamics.

Isothermal flow of two components
The entropy production in Eq.15 has two terms when two immiscible components flow at isothermal conditions.
We choose the formulation that has variables J V and J D ; volume flux and interdiffusion flux, respectively. Equation 43 gives then: where L ij ≡ L ij /T . The coefficients reflect, as above, the mechanism of flow (pressure, diffusion). Four experiments can be done to determine the four coefficients. There are only three independent coefficients. When four experiments are done, we can check the Onsager relations.

The hydraulic permeability
The (hydraulic) permeability is a main coefficient. By introducing the driving force for the volume flow from 36, we obtain With the present definition of variables the equation applies to the overall behavior of the system. A plot of J V vs.
∆p may show a threshold. This threshold has more contributions than in the single component system, as there are contributions to the pressure from surface and line energies. A threshold may be detectable at low capillary numbers.
The hydraulic or volume permeability, L pp , is found by measuring the volume flow caused by the overall pressure difference at uniform composition; The coefficient is a function of the saturation L pp = L pp (p, T, η,Ŝ w ). In the hydrodynamic regime, the coefficient can be modeled, assuming Poiseuille flow and the effective viscosity η eff = η wŜw + η nŜn [21].

The interdiffusion coefficient
The main coefficient L dd is an interdiffusion coefficient. It is defined at uniform pressure from the difference flux created by a difference in saturation; where we used Eq.22 for the driving force.

The coupling coefficients
The coupling coefficients in Eqs.46 express that a separation of components can be caused by a pressure gradient (L dp ) and that a volume flow can be promoted by a gradient in saturation (L pd ).
Consider first the determination of L dp . A pressure gradient may build as a consequence of a difference in composition [24]. The volume flux continues until a balance of forces is reached: From the force-balance across the system, we obtain: This condition can be used to find the unknown coupling coefficient, once the hydraulic permeability is known.
The remaining coupling coefficient can be found from the flux ratio, r, that has been called the reflection coefficient r, see F/orland [24]. At constant saturation, we have We are now in a position to compare L pd and L dp . The state of the system must be (approximately) the same, when the comparison is made.

Non-isothermal flow of two components
The set of equations 43 describe non-isothermal flow in porous media. The coefficients, L pp , L pd = L dp , L dd in the lower right-hand side corner of the conductivity matrix, were discussed above. The new coefficients are those related to heat transport. The coefficient L qq represent the Fourier type heat conductivity at uniform composition and pressure. The coefficients L pq and L dq are coupling coefficients.
cold water into warm reservoirs may thus lead to separation. Likewise, a pressure difference can arise from a temperature difference. This is called thermal osmosis [25].
Separation caused by a thermal driving force was observed in clay-containing soils where water was transported in clay capillaries against a pressure gradient. The coefficient, measured at constant pressure, was called the segregation potential [26]. The coefficient can be obtained from Eq. 43, setting ∆p = 0 and ∆µ c w = 0 (∆Ŝ w ) in the second line. We obtain This coefficient can also be found from steady state conditions, when the thermal gradient is balanced by a gradient in saturation (chemical potential) Determination of L pq requires knowledge of L pd .
The coupling coefficient L qp can be found by measuring the heat flux that accompanies the volume flux for constant composition and at isothermal conditions.
These effects have, to the best of our knowledge, so far not been measured for porous media. For homogeneous media they are well known.

Discussion and Conclusion
We have used the new formulation of the entropy production [1] to find constitutive equations for flow of two immiscible fluids in a porous medium under uniform or varying temperature, pressure and composition. Several of the equations are new in the context of porous media, but they follow well documented tracks in classical nonequilibrium thermodynamics [1,[22][23][24]. Experimental observations exist on single fluid flow, that give support to the theoretical description.
We have defined the pressure of the REV, and used the definition to define the pressure part of the driving force. The force obtains contributions from the surface -and, in principle also line -tensions of the system. This distinguishes the present formulation from their counterpart for homogeneous systems [1,[22][23][24]. We have seen that surface contributions can be spelled out for varying conditions, under the assumptions that the additive properties of the REV are Euler homogeneous of the first order. Doing this, we have been able to explain for instance deviations from Darcy's law, the occurrence of threshold pressures. We have defined the transport equations on the macroscopic scale, and pointed at possibilities to describe non-isothermal phenomena. As for instance Eq.27 shows, there is a multitude of scenarios that can be further investigated, and used to check the theory, cf. subsections 5. 1-5.3. The expressions open up the possibility to test the thermodynamic models in use, for their compatibility with the second law.
The basic assumption that the REV set of basis variables are Euler homogeneous functions of degree one means in its essence that one temperature, one pressure and one chemical potential per component can be defined for the REV. The properties can vary in time and space, and can therefore also be used for a transient description. We did not consider surface areas or their curvature as independent variables, but these may be included as variables. In that sense, our description be equivalent to a description using Minkovsky functionals [14]. The state of the REV is a thermodynamic state, meaning that the REV is ergodic. Some evidence already supports this idea [4], [5], originally proposed by Hansen and Ramstad [3] and Tallakstad and coworkers [16].
We have often used a specific case to illustrate relations; the non-isothermal flow of one or two immiscible fluids in a non-deformable medium. It is straight forward to include more terms in the chemical potential (e.g. gravity).
To include stress fields or other fields is more problematic.
Flow of two isothermal, immiscible fluids in a porous medium has often been described by Darcy's law, using the relative permeability concept. The seepage velocities v n and v w are related to fluxes used here by v n = J n V n and v w = J w V w . The expressions for the seepage velocities must be contained or be equivalent to the expressions given here, using the condition of invariance for the entropy production. A comparison can elucidate assumptions that are made. Hilfer and Standnes et al. [27,28] gave a set of linear relations for the seepage velocities. Their driving forces were the gradients in the single component pressures, obtained by pressure measurements in the single phases. Their description implies e.g that the composition is uniform.
We have seen through these examples how non-equilibrium thermodynamic theory gives a fundamental basis to the constitutive equations. By demanding that transport properties obey entropy production invariance and Onsager symmetry, we can also find relations between variables used so far.