Non-isothermal transport of multi-phase fluids in porous media. The entropy production

We derive the entropy production for transport of multi-phase fluids in a non-deformable, porous medium exposed to differences in pressure, temperature, and chemical potentials. Thermodynamic extensive variables on the macro-scale are obtained by integrating over a representative elementary volume (REV). Using Euler homogeneity of the first order, we obtain the Gibbs equation for the REV. From this we define the intensive variables, the temperature, pressure and chemical potentials and, using the balance equations, derive the entropy production for the REV. The entropy production defines sets of independent conjugate thermodynamic fluxes and forces in the standard way. The transport of two-phase flow of immiscible components is used to illustrate the equations.


Introduction
The aim of this article is to develop a macro-scale description of multi-phase flow in porous media in terms of non-equilibrium thermodynamics. The system consists of several fluid phases in a medium of constant porosity.
The aim is to describe transport on the scale of measurements; i.e. on the macro-scale, using properties defined on this scale, which represent the underlying structure on the micro-scale. The effort is not new; it was pioneered more than 30 years ago [1][2][3][4], and we shall build heavily on these results, in particular those of Hassanizadeh and Gray [2,3]. The aim is also still the original one; to obtain a systematic description, which can avoid arbitrariness and capture the essential properties of multi-component multi-phase flow-systems. Not only bulk energies need be taken into account to achieve this. Also the excess surface-and line-energies must be considered.
But, unlike what has been done before, we shall seek to reduce drastically the number of variables needed for the description, allowing us still to make use of the systematic theory of non-equilibrium thermodynamics. While the entropy production in the porous medium so far has been written as a combination of contributions from each phase, interface and contact line, we shall write the property for a more limited set of macro-scale variables. This will enable us to describe experiments and connect variables at this scale.
The theory of non-equilibrium thermodynamics was set up by Onsager [5,6] and further developed for homogeneous systems during the middle of the last century [7]. It was the favored thermodynamic basis of Hassanizadeh and Gray for their description of porous media. These authors [2,3] discussed also other approaches, e.g the theory of mixtures in macroscopic continuum mechanics, cf. [1,4].
The theory of classical non-equilibrium thermodynamics has been extended to deal with a particular case of flow in heterogeneous systems, namely transport along [8] and perpendicular [9] to layered interfaces. A description of heterogeneous systems on the macro-scale has not been given, however. Transport in porous media take place, not only under pressure gradients. Temperature gradients will frequently follow from transport of mass, for instance in heterogeneous catalysis [10], in polymer electrolyte fuel cells, in batteries [9,11], or in capillaries in frozen soils during frost heave [12]. The number of this type of phenomena is enormous. We have chosen to consider only the vectorial driving forces related to changes in pressure, chemical composition and temperature, staying away for the time being from deformations, chemical reactions, or forces leading to stress [13]. The multi-phase flow problem is thus in focus.
The development of a general thermodynamic basis for multi-phase flow [2,3] started by introduction of thermodynamic properties for each component in each phase, interface and three-phase contact line. A representative volume element (REV) was introduced, consisting of bulk phases, interfaces and three-phase contact lines. Balance equations were formulated for each phase in the REV, and the total REV entropy production was the sum of the separate contributions from each phase.
In a recent work by Hansen et al. [14], it was recognized that the description of the motion of the fluids at the coarse-grained level could be described by extensive variables. The properties of Euler homogeneous functions could then be used to create relations between the flow rates at this level of description. This work, however, did not address the coarse-graining problem itself. We shall take advantage of Euler homogeneity also here and use it in the coarse-graining process described above.
Like Hassanizadeh and Gray [2,3] we use the entropy production as the governing property. But rather than dealing with the total entropy production as a sum of numerous parts, we shall seek to define the total entropy production directly from a basis set of a few coarse-grained variables. Central in the work is the REV, essentially the same as defined before. With the REV as a starting point we will derive the entropy production in this part of the work. The independent and measurable variables are defined here. The independent constitutive equations that follow from the entropy production will be given in a next part of the work. Doing this, we aim to contribute towards solving the scaling problem; i.e. how a macro-level description can be obtained consistent with the micro-level one. The Figure 1: Schematic illustration of a heterogeneous system (white box, length L) exposed to a difference in temperature, ∆T , pressure, ∆p or chemical potential ∆µ i . The system is isolated in the y, z−directions. Net flows take place in the x−direction. Three representative elementary volumes, REVs (magenta squares, length l) are indicated. The REVs may overlap. Each is represented by a set of variables (p, T, µ i ) which defines a state (blue dot). Such states can be defined anywhere on the x-axis.
overall aim is to describe transport due to thermal, chemical, mechanical and gravitational forces in a systematic, course-grained manner that is sufficiently simple for practical use.

System
Consider a heterogeneous system as illustrated by the (white) box in Fig. 1. The system is a porous medium of fixed porosity filled with several immiscible fluids. There is net transport in one direction only, the x-direction. On the scale of measurement, the system is without structure. By zooming in, we see the pore scale. A collection of pores with two fluids is schematically shown in Figure 2.
A temperature, pressure or chemical potential difference is applied between the inlet and the outlet, and these differences can be measured. The pressure difference ∆p between the outlet and the inlet was defined for steady state conditions by Tallakstad et al. [15], as the time average of the fluctuating difference ∆p(t): Here t is the time. Subscript 'b' denotes beginning and 'e' denotes the end of the measurement. We adopt similar definitions for ∆T and ∆µ i . It is possible, through application of separate inlet channels, to control the flow into and out of the system and find the flow of each component, to define the flow situation in Fig. 1. In the presence of two immiscible phases, it is only possible to define the pressure difference between the inlet and the outlet for the phases, ∆p w and ∆p n , if there is continuity in the respective phases. We will repeatedly use two-phase flow as example, where w indicates the most wetting and n the least wetting phase. We refer to them simply as the wetting and the non-wetting phase. In most of the paper we consider a multi-phase fluid. In the system pictured in Fig. 2, there is flow within the REV in the direction of the pore, which we will also call a link. This is not necessarily the direction given by the overall pressure gradient. The flow on the macro-scale, however, is always in the direction of the pressure gradient. Net flow in other directions are zero due to isolation of the system in these directions. By flow on the macro-scale, we mean flow in the direction of the overall pressure gradient along the x-coordinate in Fig.1. The value of this average flow is of interest.
The representative volume element, REV, is constructed from a collection of pores like those contained in the red square in Fig.2. In Fig.1, three REVs are indicated (magenta structured squares). In a homogeneous system, statistical mechanical distributions of molecular properties lead to the macroscopic properties of a volume element.
In a heterogeneous system like here, the statistical distributions are over the states within the REV. The collection of pores in the REV, cf. Fig. 2, should be of a size that is large enough to provide meaningful values for the extensive variables, and therefore well defined intensive variables (see below, Eqs. 18 and 19), cf. section 3.2 below. Thermodynamic relations can be written for each REV.
State variables characterize the REV. They are represented by the (blue) dots in Fig. 1. The size of the REV depends on its composition and other conditions. Typically, the extension of a REV, l, is large compared to the pore size of the medium, and small compared to the full system length L. This construction of a REV is similar to the procedure followed in smoothed particle hydrodynamics [16], cf. the discussion at the end of the work.
The REVs so constructed, can be used to make a path of states, over which we can integrate across the system.
Each REV in the series of states, is characterized by variables T, p, µ i , as indicated by the blue dots in Fig.1. Vice versa, each point in a porous medium can be seen as a center in a REV. The states are difficult to access directly, but can be accessed via systems in equilibrium with the states, as is normal in thermodynamics. This is discussed at the end of the work. We proceed to define the REV-variables.
3 Properties of the REV

Porosity and saturation
Consider a solid matrix of constant porosity φ. We are dealing with a class of systems that are homogeneous in the sense that the typical pore diameter and pore surface area are, on the average, the same everywhere. The pores are filled with a mixture of m − 1 fluid phases; the solid matrix is labeled m. Properties will depend on the time, but this will not be indicated explicitly in the equations.
The system is filled with m phases. In a simple case, the phases are immiscible, and chemical constituents are synonymous with a phase. The state of the REV can be characterized by the volumes of the fluid phases V α,REV , α = 1, .., m − 1 and of the solid medium V m,REV . The total volume of the pores is while the volume of the REV is Superscript REV is used to indicate a property of the REV. The last term is the sum of the excess volumes of the three-phase contact lines. While the excess volume of the surfaces is zero by definition, this is not the case for the three-phase contact lines. The reason is that the dividing surfaces may cross each other at three lines which have a slightly different location. The corresponding excess volume is in general very small, and will from now on be neglected. This gives the simpler expression All these volumes can be measured.
The porosity, φ, and the saturation,Ŝ, are given by The porosity and the saturation are intensive variables. They do not depend on the size of the REV. They have therefore no superscript. It follows from these definitions that In addition to the volumes of the different bulk phases (they are fluids or solids) m ≥ α ≥ 1, there are interfacial areas, Ω, between each two phases in the REV: Ω αβ,REV , m ≥ α > β ≥ 1. The total surface area of the pores is measurable. It can be split between various contributions When the surface is not completely wetted, we can estimate the surface area between the solid m and the fluid phase α, from the total pore area available and the saturation of the component.
This estimate is not correct for strongly wetting components or dispersions. In those cases, films can be formed at the walls, and Ω mα,REV is not proportional toŜ α . In the class of systems we consider, all fluids touch the wall, and there are no films of one fluid between the wall and another fluid.

Thermodynamic properties of the REV
We proceed to define the thermodynamic properties of the REV within the volume V REV described above. In addition to the volume, there are other additive variables. They are the masses, the energy and the entropy. We label the components (the chemical constituents) using italic subscripts. There are in total n components distributed over the phases, surfaces and contact lines. The mass of component i, M REV i , in the REV is the sum of bulk masses, , and the excess line masses, There is some freedom in how we allocate the mass to the various phases and interfaces [9,17]. We are e.g. free to choose a dividing surface such that one M αβ,REV i equals zero. A zero excess mass will simplify the description, but will introduce a reference. The dividing surface with zero M αβ,REV i is the equimolar surface of component i.
The total mass of a component in the REV is, however, independent of the location of the dividing surfaces. From the masses, we compute the various mass densities where ρ i and ρ α i have dimension kg.m −3 , ρ αβ i has dimension kg.m −2 and ρ αβδ i has dimension kg.m −1 .
All densities are for the REV. If we increase the size of the REV, by for instance doubling its size, V REV , M REV i and other extensive variables will all double. They will double, by doubling all contributions to these quantities.
But this is not the case for the density ρ i or the other densities. They remain the same, independent of the size of the REV. This is true also for the densities of the bulk phases, surfaces and contact lines. Superscript REV is therefore not used for the densities.
Within one REV there are natural fluctuations in the densities. But the densities make it possible to give a description on the macro-scale independent of the precise size of the REV. The densities will thus be used in our description on the macro-scale. The density ρ α i may vary somewhat in V α . We can then find M α i as the integral of ρ α i over V α . Equation 10 then gives the volume-averaged densities.
The internal energy of the REV, U REV , is the sum of bulk internal energies, U α,REV , m ≥ α ≥ 1, the excess interfacial internal energies, U αβ,REV , m ≥ α > β ≥ 1, and the excess line internal energies, U αβδ,REV , m ≥ α > β > δ ≥ 1: The summation is taken over all phases, interfaces and contact lines (if non-negligible). We shall see in a subsequent paper how these contributions may give specific contributions to the driving force. The internal energy densities are defined by The entropy in the REV, S REV , is the sum of bulk entropies, S α,REV , m ≥ α ≥ 1, the excess entropies, S αβ,REV , m ≥ α > β ≥ 1, and the excess line entropies, S αβδ,REV , m ≥ α > β > δ ≥ 1: The entropy densities are defined by and To summarize this section; we have defined a basis set of variables for a class of systems, where these variables are additive in the manner shown. From the set of REV variables we obtain the densities, u, s or ρ i to describe the heterogeneous system on the macro-scale. A series of the REVs of this type, is needed for integration across the system, see section 5.
The surface areas and the contact line lengths are not independent variables in this representation of the REV.
These variables are, however, included indirectly through the assumption that the basic variables of the REV are additive. This means that a REV of a double size has double the energy, entropy, and mass, but also double the surface areas of various types and double the line lengths.

REV size considerations
As an illustration of the REV construction, consider the internal energy of two isothermal, immiscible and incompressible fluids (water and decane) flowing in a Hele-Shaw cell composed of silicone glass beads. The relevant properties of the fluids can be found in Table 1. The porous medium is a hexagonal network of 3600 links, as illustrated in Figure 3. The network is periodic in the longitudinal and the transverse directions and a pressure difference of 1.8 × 10 4 Pa drives the flow in the longitudinal direction. The overall saturation of the water is 0.4.
The network flows were simulated using the method of Aker et al. [19], see Ref. [20] for details.
The internal energy of the REV is, according to Section 3.2, a sum over the two fluid bulk contributions and three interface contributions, where, u i is the internal energy density of phase i and u ij is the excess internal energy per interfacial area between phase i and phase j. We assume u i and u ij to be constant. For simplicity, u ij is approximated by interfacial tension, denoted γ ij . The internal energy of the porous matrix is constant in this example and is therefore set to zero.
Candidate REVs are of different sizes, see Table 2. The 5.4 mm × 6 mm (green), and 10.4 mm × 12 mm (blue) candidate REVs are shown in Figure 3. For all candidate REVs, U REV is calculated according to Eq. 16 at each time step. Since the measured saturations and interfacial areas are fluctuating in time, so is the internal energy. A time-step weighted histogram of the internal energy presents the probability distribution.
The probability distributions of U REV are shown in Figure 4 for the 5.2 mm × 6 mm (green) and 10.4 mm × 12 mm (blue) candidate REVs. In both plots, the vertical lines represent the internal energy the REV would have if it were occupied by one of the fluids alone. We denote the difference in internal energy between these two single-phase states by ∆U REV .
The mean value of the U REV for all candidate REVs are given in Table 2, along with mean density u = U REV /V REV and the standard deviation of U REV divided by ∆U REV . The latter quantity is a measure of the relative size of the   This example indicates that it makes sense to characterize the internal energy of a porous medium in terms of an internal energy density as defined by Eq. 11, given that the size of the REV is appropriately large. About 100 links seem to be enough in this case. This will vary with the type of porous medium, cf. the 2D square network model of Savani et al. [24].

Homogeneity on the macro-scale
Before we address any transport problems, consider again the system pictured in Fig.1 (the white box). All REVs have variables and densities as explained above. By integrating to a somewhat larger volume V , using the densities defined, we obtain the set of basis variables, (U, S, M i ), in V . The internal energy U of the system is an Euler homogeneous function of first order in S, V, M i : The internal energy U, volume V , entropy S, and component mass M i , obey therefore the Gibbs equation; No special notation is used here to indicate that U, S, V, M i are properties on the macro-scale. Given the heterogeneous nature on the micro-scale, the internal energy has contributions from all parts of the volume V , including from the excess surface and line energies. By writing Eq.17 we find that the normal thermodynamic relations apply for the heterogeneous system at equilibrium, additive properties U, S, V, M i , obtained from sums of the bulk-, excess surface-and excess line-contributions.
We can then move one more step and use Gibbs equation to define the temperature, the pressure and chemical potentials on the macro-scale as partial derivatives of U : The temperature, pressure and chemical potentials on the macro-scale are, with these formulas, defined as partial derivatives of the internal energy. This is normal in thermodynamics, but the meaning is now extended. In a normal homogeneous, isotropic system at equilibrium, the temperature, pressure and chemical equilibrium refer to a homogeneous volume element. The temperature of the REV is a temperature representing all phases, interfaces and lines combined, and the chemical potential of i is similarly obtained from the internal energy of all phases.
Therefore, there are only one T, p and µ i for the REV. The state can be represented by the (blue) dots in Fig.1.
On the single pore level, the pressure and temperature in the REV will have a distribution. In the two immisciblefluid-example the pressure, for instance, will vary between a wetting and a non-wetting phase because of the capillary pressure. One may also envision that small phase changes in one component (e.g. water) leads to temperature variations due to condensation or evaporation. Variations in temperature will follow changes in composition.
The intensive properties are not averages of the corresponding entities on the pore-scale over the REV. This was pointed out already by Gray and Hassanizadeh [3]. The definitions are derived from the total internal energy only, and this makes them uniquely defined. It is interesting that the intensive variables do not depend on how we split the energy into into bulk and surface terms inside the REV.
By substituting Eq.19 into Eq.18 we obtain the Gibbs equation for a change in total internal energy on the macro- As a consequence of the condition of homogeneity of the first order, we also have The partial derivatives T, p and µ i are homogeneous functions of the zeroth order. This implies that Choosing λ = 1/V it follows that The temperature therefore depends only on the subset of variables s ≡ S/V, ρ i ≡ M i /V and not on the complete set of variables S, V, M i . The same is true for the pressure, p, and the chemical potentials, µ i . This implies that T, p and µ i are not independent. We proceed to repeat the standard derivation of the Gibbs-Duhem equation which makes their interdependency explicit.
The Gibbs equation on the macro-scale in terms of the densities follows using Eqs. 20 and 21 which can alternatively be written as The Euler equation implies By differentiating Eq. (26) and subtracting the Gibbs equation (24), we obtain in the usual way the Gibbs-Duhem equation: This equation makes it possible to calculate p as a function of T and µ i and shows how these quantities depend on one another.
We have now described the heterogeneous porous medium by a limited set of coarse-grained thermodynamic variables. These average variables and their corresponding temperature, pressure and chemical potentials, describe a coarse-grained homogeneous mixture with variables which reflect the properties of the class of porous media. In standard equilibrium thermodynamics, Gibbs' equation applies to a homogeneous phase. We have extended this use to be applicable for heterogeneous systems at the macro-scale. Whether or not the chosen procedure is viable, remains to be tested. We refer to the Section 7 of this paper.

Entropy production in porous media
Gradients in mass-and energy densities produce changes in the variables on the macro-scale. These lead to transport of heat and mass. Our aim is to find the equations that govern this transport across the REV. We therefore expose the system to driving forces and return to Fig. 1.
The balance equations for masses and internal energy of a REV are The transport on this scale is in the x−direction only. The mass fluxes, J i , and the flux of internal energy, J u , are all macro-scale fluxes. The internal energy flux is the sum of the measurable (or sensible) heat flux, J q and the partial specific enthalpy (latent heat), H i (in J.kg −1 ) times the component fluxes, J i , see [3,7,9] for further explanations. Component m (the porous medium) is not moving and is the convenient frame of reference for the fluxes.
The entropy balance on the macro-scale is Here J s is the entropy flux, and σ is the entropy production which is positive definite, σ ≥ 0 (the second law of thermodynamics). We can now derive the expression for σ in the standard way [7,9], by combining the balance equations with Gibbs' equation. The entropy production is the sum of all contributions within the REV.
In the derivations, we assume that the Gibbs equation is valid for the REV also when transport takes place. Droplets can form at high flow rates, while ganglia may occur at low rates. We have seen above that there is a minimum size of the REV, for which the Gibbs equation can be written. When we assume that the Gibbs equation applies, we implicitly assume that there exists a uniquely defined state. The existence of such a state was first postulated by Hansen and Ramstad [25]. Experimental evidence for the assumption was found by Erpelding [26]. giving We can now introduce the balance equations for mass and energy into this equation, see [9] for details. By comparing the result with the entropy balance, Eq.29, we identify first the entropy flux, J s , The entropy flux is composed of the sensible heat flux over the temperature plus the sum of the specific entropies carried by the components. The form of the entropy production, σ, depends on our choice of the energy flux, J u or J q . The choice of form is normally motivated by practical wishes; what is measurable or computable. We have

The chemical potential at constant temperature
The derivative of the chemical potential at constant temperature is needed in the driving forces in the second line for σ in Eq.32. For convenience we repeat its relation to the full chemical potential [7]. The differential of the full chemical potential is: where S i , V i and (∂µ i /∂M j ) p,T,Mi are partial specific quantities. The partial specific entropy and volume are equal to: and the last term of Eq. 32 is denoted by By reshuffling, we have the quantity of interest as the differential of the full chemical potential plus an entropic term; The differential of the chemical potential at constant temperature is With equilibrium in the gravitational field, the pressure gradient is dp/dx = −ρg, where ρ is the total mass density and g is the acceleration of free fall [27]. The well known separation of components in the gravitational field is where W i is the molar mass (in kg.mol −1 ),Ŝ i the saturation, and y i the activity coefficient of component i.
The gas constant, R, has dimension J.K −1 .mol −1 . The gradient of the mole fraction of methane and decane in the geothermal gradient of the fractured carbonaceous Ekofisk oil field, was estimated to 5x10 −4 m −1 [28], in qualitative agreement with observations. We replace dµ i,T below using these expressions.
It follows from Euler homogeneity that the chemical potentials in a (quasi-homogeneous) mixture are related by 0 = SdT − V dp + n j=1 ρ j dµ j , which is Gibbs-Duhem's equation. By introducing 36 into this equation we obtain an equivalent expression, to be used below: 6 Transport of heat and two-phase fluids Consider again the case of two immiscible fluids, one more wetting (w) and one more non-wetting (n). The entropy production in Eq.32 gives, The solid matrix is the frame of reference for transport, J r = 0 and does not contribute to the entropy production.
The volume flux is frequently measured, and we wish to introduce this as new variable Here J V has dimension (m 3 .m −2 .s −1 = m.s −1 ), and the partial specific volumes have dimension m 3 .kg −1 . The chemical potential of the solid matrix may not vary much if the composition of the solid is constant across the system. We assume that this is the case (dµ c m ≈ 0), and use Eq.39 to obtain expression for dµ i,T into Eq.40, and find the practical expression: In the last line, the difference velocity J D is This velocity describes the relative movement of the two components within the porous matrix. In other words, it describes the ability of the medium to separate components. The main driving force for separation is the chemical driving force, related to the gradient of the saturation. The equation implies that also temperature and pressure gradients may play a role for the separation.
The entropy production has again three terms, one for each independent driving force. With a single fluid, the number of terms are two. The force conjugate to the heat flux is again the gradient of the inverse temperature. The entropy production, in the form we can obtain, Eqs.40 or 43, dictates the constitutive equations of the system.
The volume flows used by Hansen et al [14] are related to ours by J n V n = Q n /A p and J w V w = Q w /A p , where A p ≡ V p / is the cross-sectional area of the pores in the volume normal to the flow direction, and is the length of the volume in the flow direction. In [14] the total volume flow was assumed to be an Euler homogeneous functions of the first order of the fractional areas, A n ≡ V n / and A w ≡ V w / for the two components, respectively. The seepage velocities are v n = Q n /A n and v w = Q w /A w . By introducing these variables in the expression for the entropy production, we obtain two alternative forms for the entropy production, Eq. 43: In the second identity we used Gibbs-Duhem, Eq. 42. If we define the entropy production can be written as Whether p w and p n have any concrete meaning is not clear. The difference of their gradients is where we again used Gibbs-Duhem in the last identity. As one often identifies the left hand side of this equation with the gradient of the capillary pressure, this would relate this force to the chemical force. We will come back to this issue later. Figure 5: A one-dimensional heterogeneous system cut into slices. Each cut is brought in equilibrium with a homogeneous mixture at the same temperature and pressure as the REV.

A path of sister systems
As pointed out above, through the construction of the REV we were able to create a continuous path through the system, defined by the thermodynamic variables of the REVs. The path was illustrated by a sequence of dots in Fig. 1. Such a path must exist, to make integration possible. Also continuum mixture theory hypothesizes such a path [4]: Hilfer introduced a series of mixture states, to define an integration path across the porous system, see e.g. [4].
The path created in Section 2 is sufficient as a path integration across the medium. The access to and measurement of properties in the REVs is another issue. It is difficult, if not impossible, to measure in situ as stated upfront. The measurement probe has a minimum extension (of some mm), and the measurement will represent an average over the surface of the probe. For a phase with constant density, the average is well defined and measurable. A link between the state of the REV and a state where measurements are possible, is therefore needed. We call the state that provides this link a sister state.
Consider again the path of REVs in the direction of transport. To create the link between the REV and its sister state, consider the system divided into slices, see Figure 5. The slice (the sister system) contains homogeneous phases in equilibrium with the REV at the chosen location.
We hypothesize that we can find such sister states; in the form of a multi-component mixture with temperature, pressure and composition such that equilibrium can be obtained with for the REV-variables at any slice position.
The variables of the sister state can then be measured the normal way. The chemical potential of a component in the sister state can, for instance, be found by introducing a vapor phase above this state and measure the partial vapor pressure. We postulate thus that a sister state can be found, that obey the conditions Here i=1,...n are the components in the REV, and superscript s denotes the sister state. With the sister states available, we obtain an experimental handle on the variables of the porous medium. The hypothesis must be checked, of course.
The series of sister states have the same boundary conditions as the REV-states, by construction, and the overall driving forces will be the same. Between the end states, we envision the non-equilibrium system as a staircase.
Each step in the stair made up of REV is in equilibrium with a step of the sister-state-stair. Unlike the states inside the porous medium, the sister states are accessible for measurements, or determination of T, p and µ i . Transport can then be described as taking place along a sequence of the sister states.

Discussion
We have seen how we can define a representative elementary volume, called the REV, and conditions such that standard non-equilibrium thermodynamics for homogeneous systems [7] can be extended to describe flow in heterogeneous systems like porous media. The procedure adds a step to earlier elaborations in the field [2,3]. We have been able to derive the total entropy production directly for the heterogeneous system, and have thus obtained a new basis for constitutive equations of the system on the macro-scale. Earlier, the equation sets were written for the separate parts, or for systems that were homogeneous on the molecular scale [29].
Equations 40 and 43 have been established since long in the theory of transport in polymer membranes, see e.g. [29]. The step to include porous media is inspired by this. The way of dealing with lack of knowledge of variables inside the system was for instance used in polymer membrane transport long ago, see [30,31]. The procedure, to introduce a series of equilibrium states, each state in equilibrium with the membrane at some location between the external boundaries, was first used by Scatchard [30], and later experimentally verified [31,32].
The clue to obtain the present results was to use Euler homogeneity on the macro-scale. This allowed us to set up the Gibbs equation with well defined variables, in spite of local variations in, for instance, the pressure inside the REV. Together with the balance equations, this gave the entropy production in a standard way. Some support can be found in the literature. Prigogine and Mazur [33] investigated a mixture of two fluids using non-equilibrium thermodynamics. Their system consisted of superfluid -and normal helium. Two different pressures were defined.
For the case they discussed, the interaction between the two fluids was small, meaning that one phase flowed as if the other one (aside from a small frictional force) was not there. The situation here is similar, as we may have different liquid pressures inside the REV. But the interaction between the two immiscible components in our porous medium is large, not negligible as in the helium case.
A central issue in the extension is the construction of the REV. This construction is subject to conservation of the important properties of the system. For instance, the REV must conserve the total porosity of the system. Such constructions are known from other up-scaling techniques. In mesoscale simulations, i.e. upscaling from properties at the nanometer scale to micrometer scale, techniques known as "Dissipative Particle Dynamcs" stemming from hydrodynamics [34] and "Smoothed Particle Hydrodynamics", stemming from astrophysics [16] have been used.
Inspired by the idea behind smoothed particle hydrodynamics, we can also define a normalized weight function W (r), such that a microscopic variable a(r) may be represented by its average, defined as For example, if a(r ) is the local void fraction in a porous material as determined from samples of the material, a(r) is the average porosity of the medium. The average is assigned to the point r and varies smoothly in space.
The average porosity a(r) would then be suitable for e.g. a reservoir simulation at the macro-scale.
In general, the system is subject to external forces and its properties are non-uniform. The choice of W (r) is therefore crucial in that it defines the extent of the coarse-graining and the profile of the weighting. The illustration in Fig. 1 alludes to a weight function that is constant inside a cubic box and zero outside, but other choices are possible. Popular choices used in mesoscale simulations are the Gaussian and spline functions (See [16] for details).
A convenient feature of the coarse-graining is that the average of a gradient of a property a is equal to the gradient of the average. ∇a(r) = ∇a(r) Similarly the average of a divergence of a flux is equal to the divergence of the average. This implies that balance equations, which usually contain the divergence of a flux, remain valid after averaging.
Time averages can also be introduced along the lines sketched above. Also the average of a time derivative of a becomes equal to the time derivative of the average. Again, the balance equations remain valid after averaging. (In this paper we have not used time averages.) The total amount of a in the volume of the REV is The time scale relevant to porous media transport are usually large (minutes, hours); and much larger than times relevant for the molecular scale. Partial derivatives were used in the balance equations and the entropy production.
This signifies that properties can change not only along the coordinate axis, but also on the time scale. In the present formulation any change brought about in the REV must retain the validity of the Gibbs equation (i.e. Euler homogeneity). As long as that is true, we may use the equations, also for transient phenomena.
The outcome of the derivations enables us to deal with a wide range of non-isothermal phenomena in a systematic manner, from frost heave to heterogeneous catalysis, or multi-phase flow in porous media. We will elaborate on what this means in the next part of this work. In particular, we shall give more details on the meaning of the additive variables and the REV pressure. We will return to the meaning of the REV variables and how they will contribute to the driving forces of transport.

Concluding remarks
We have derived the entropy production for transport of heat and mass in porous media. The derivations have followed standard non-equilibrium thermodynamics for heterogeneous systems [9]. The only, but essential, difference has been the fact that we write all these equations for a porous medium on the macro-scale for the REV of a minimum size using its total entropy, energy and mass, while these equations often are only written for the separate contributions. Broadly speaking, we have been zooming out our view on the porous medium to first define some states that we take as thermodynamic states because they obey Euler homogeneity. The states are those illustrated by the dots in Figure 1. In order to define these states by experiments, we constructed the sister states of Fig.5.
The advantage of the present formulations is this; it is now possible to formulate the transport problem on the scale of a flow experiment in accordance with the second law of thermodynamics, with far less variables. This also opens up the possibility to test the thermodynamic models for use along with and compatibility with this law.
Such tests can be explicitly formulated once we have the constitutive equations. In the next part of this work, such details will be given.
We have been dealing as an example with the specific case of immiscible, non-isothermal two-phase flow in a non-deformable medium. In deformable systems, grains may move by translation, compaction and/or rotation. We expect that this situation also can be treated by the procedure given above.