Fluctuation-dissipation theorems for flow in porous media

A thermodynamic description of nano-porous media must handle the sizeand shape-dependence of the media properties. Such dependencies are typically due to the presence of immiscible phases, contact areas and contact lines. We propose a way to obtain average densities suitable for integration on the course grained scale, applying Hill’s thermodynamics for small systems to the subsystems. we argue that the average densities of the porous medium, when defined in a proper way, obey the Gibbs equation. All contributions are additive or weakly coupled. From the Gibbs equation and the balance equations, we derive the entropy production in the standard way, for transport of multi-phase fluids in a non-deformable, porous medium exposed to differences in boundary pressures, temperatures, and chemical potentials. Linear relations between thermodynamic fluxes and forces follow for the control volume. Fluctuationdissipation theorems are formulated for the first time, for the fluctuating contributions to fluxes in the porous medium. These give an added possibility for determination of porous media permeabilities. Practical possibilities are further discussed.


Introduction
Porous media represent a vast and important class of systems; present for instance in biology, geology or in technological applications. Our interest is to describe transport in porous media. Examples are transport of nanoparticles in cancerous tissue, migrations in ground water reservoirs, or heterogeneous catalysis. A description of porous media transport must reflect on the macro-scale what goes on on the smaller scales. The choice of variables on the continuum scale is therefore central. For porous media, it has been difficult to find good bottom-up descriptions of the complex, heterogeneous structure of the medium; as it contains a mixture of solids and often immiscible fluids, wetting or non-wetting, microporous and nanoporous. Methods that describe fluids, confined in nanopores, are also scarce in spite of the documented large impact of confinement [12].
Porous media may be exposed to several fields. Consider the case of a porous reservoir rock. On the macro-scale, we may have a variation in pressure in addition to the gravitational field. A gradient in composition may also be present. In the fractured carbonaceous Ekofisk oil field a geothermal gradient was found to be relevant [13]. The interplay of several driving forces is then central. A suitable thermodynamic theory of transport must be able to describe such interplay and the energy dissipation involved (the entropy production). The fluctuations in each flux will be yet a characteristic signature of the medium.
Thermodynamic theories of transport in homogeneous and heterogeneous fluids and solids have been studied over most of the last century and are well documented in nonequilibrium thermodynamics, see e.g. [1][2][3]. Also hydrodynamics fluctuations are well studied for bulk fluids [4]. Transport processes in porous media occur on a wider range of time scales. Winkler et al. [29] observed athermal relaxation times varying from 10 −3 to 0.1 seconds for viscous flow in a network. By including coalescence dynamics and interface phenomena and pore scale dynamics, the scale of Winkler et al. need be extended at both ends, down to 10 −5 and up to 1-10 seconds [32]. Fluid configuration shifts take even longer, up to minutes [31]. Phenomena that are well separated on the time scale can be regarded as not coupled, and dealt with separately. In order for various phenomena to interact or couple, they need to relax on (approximately) the same time scale.
Several attempts have been made to describe multiphase flow of this kind [5][6][7][8][9], all in essence aiming at a formulation in terms of average properties, that can describe measurements or simulations. Important while the volume of the REV is V REV ≡ V r,REV + V p,REV + V nwr,REV (2) Superscript REV is used to indicate a property of the REV. The last term is the excess volume 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. The volume of the REV is a so-called additive variable. The porosity, φ, and the saturations, S n and S w , are given by The porosity and the saturations are intensive variables. In our choice of REV, 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 bulk phases, there are interfacial areas between the phases in the REV: Ω nw,REV , Ω nr,REV , Ω wr,REV . The total surface area of the solid phase is given by The total length of the three phase contact lines is given by Λ nwr,REV . The surface areas and the contact line lengths are additive. This means that a REV of a double size has double the surface areas of various types and double the line length. Two REVs have properties with double the value of one. So far our definitions are standard.

Thermodynamic properties of the REV
We proceed to define the thermodynamic properties of a REV of volume V REV following [14]. The proposal to use Minkowski functionals is similar [22,23]. We will take the REV large in the Hill sense. This implies that on the REV level, the additive variables of a REV with double the size are twice the value in the single REV. No dependence of the densities on the size of the REV. Additivity inside the REV is explained below.

Additive variables
In addition to the volume, there are other additive variables, most prominently the masses. The mass of components n, w, r in the REV is the sum of bulk, excess interfacial, and excess line masses We have chosen to use the equimolar (or equimass) surface of the rock as the fluid-rock dividing surface and the wetting fluid dividing surface as the fluid-fluid dividing surface. Furthermore, we positioned the contact lines such that M nwr,REV r = 0. The total mass of each component in the REV is independent of our choice of the location of the dividing surfaces and contact line. When we use the various mass densities Eqs.6-8 become ρ n V REV = ρ n n V n,REV + ρ nr n Ω nr,REV + ρ nw n Ω nw,REV + ρ nwr n Λ nwr,REV The mass densities have dimensions kg.m −3 for the bulk phases, kg.m −2 for the surfaces and kg.m −1 for the contact line. The masses of the components are, like the volume, additive variables.
All densities refer to the REV. The REV is chosen such that if we increase the size of the REV, by for instance doubling its size, all variables with a superscript REV will double. (The variables are additive.) But this is not the case for densities. They remain the same, independent of the size of the REV. This is in agreement with the early REV definitions [5]. All the densities of the bulk phases, surfaces and contact lines remain the same. Superscript REV is therefore not used for the densities.
The entropy of the REV is the sum of bulk, excess interfacial, excess line entropies is: The sum of bulk, excess interfacial, excess line entropies is the material contribution, S mat , to the entropy. There are many ways to distribute the phases, and the origin of the configurational contribution comes from the many configurations that have the same component volumes, surface areas and contact line lengths in the REV. The integral of the logarithm of the corresponding probability distribution times k B gives the configurational contribution to the entropy. Both the material and the configurational contributions are additive. By introducing the entropy densities, the equation becomes The dimensions are J.K −1 .m −3 for the bulk and configurational terms, J.K −1 .m −2 for the excess interfacial entropy density, and J.K −1 .m −1 for the excess line entropy density. For all variables relations similar to Eqs.12 and 13 apply. The internal energy of the REV is the sum of bulk, excess interfacial, excess line internal energies and a configurational contribution: With the internal energy densities this equation becomes uV REV = u mat + u conf V REV = u n V n,REV + u w V w,REV + u r V r,REV + u nr Ω nr,REV +u wr Ω wr,REV + u nw Ω nw,REV + u nwr Λ nwr,REV + u conf V REV Their dimensions are J.m −3 for the bulk and configurational terms, J.m −2 for the excess interfacial energy density and J.m −1 for the excess line internal energy density. The above variables (u, s, ρ n , ρ w ) constitute the basis set of thermodynamic variables that will enter the Gibbs equation for the REV.

Dealing with smallness: The grand potential and the replica energy
The REV in porous media are often open to the surroundings. The grand potential is then relevant for a description. This potential concerns an open system of fixed volume with controlled temperature and chemical potentials. The surrounding reservoir is used for control. The medium can have nano-sized pores, and the pore fluids are then no longer like bulk fluids. We shall deal with change in thermodynamic variables using the so-called Small System Method, based on Hill's original idea [30]. In this method a small representative elementary volume, V REV , is put into a bigger box. The bigger box controls the temperature and the chemical potentials. The statistical mechanical description is given by the grand canonical ensemble (GC). The grand potential of the REV is equal to minus the logarithm of the grand canonical partition function times k B T .
The partition function can be used to support the choice of variables with additive contributions (see the next section). The grand potential X GC,REV has material as well as configurational contributions. We have, similar to Eq.12 for the entropy, the grand potential X GC,REV = X mat + X conf = X n,GC,REV + X w,GC,REV + X r,GC,REV + X nr,GC,REV +X wr,GC,REV + X nw,GC,REV + X nwr,GC,REV + X conf By introducing the densities of the grand potential, this equation becomes The densities have dimension J.m −3 for a bulk and the configurational, J.m −2 for an excess interfacial and J.m −1 for the excess line grand potential densities. The grand potential densities do not depend on the size of the REV. We therefore dropped the superscript REV for the densities.
From the method of Hill [17] we define variables particular for subsystems in the REV which are small in his sense. These are systems which have considerable surface or line energies. The variables that describe the effect of smallness, or the deviation in the thermodynamical variables from the large system limit are the hat-variables in the GC ensemble. The hat-pressure, -surface tension and line tensions are given by [17,18]: These quantities depend on the size of the pores when we approach the nano-scale, and are therefore unequal to the corresponding variable without a hat. The art is to choose the REV to be so large that p, p mat and p conf do not depend on V REV . This must apply also for the nanoporous medium. For the REV, the material and the configurational pressures this implies that When we substitute Eqs.18 and 6 into Eq.17 we obtain for the REV grand potential We show below why the grand potential X GC,REV = −pV REV and as a consequence p have configurational contributions, in the case of system control by T, V REV , µ n , µ w , µ r .

Some REV-size considerations
Equations like Eq.15 can be used to estimate the relative size of bulk-, surface-and line-contributions. When we divide by V REV , we obtain The bulk contributions are proportional to V n,REV /V REV , V w,REV /V REV , V r,REV /V REV respectively, and are therefore generally of the same order of magnitude as u. The configurational contribution can be compared to an energy of mixing, and may have the same order of magnitude as u. The surface contributions are proportional to Ω nr,REV /V REV , Ω wr,REV /V REV , Ω nw,REV /V REV respectively, and therefore in general smaller than u. The line contribution is proportional to Λ nwr,REV /V REV , and is therefore in general small compared to the surface contributions and much smaller than bulk contributions. These statements about relative sizes depend on the magnitude of the excess densities. The excess mass of a surface active material can be large. Similarly, the surface tension times the area can give a sizable contribution to the internal energy.
The REV must be large compared to the single pore volume, to provide representative average values. It should not be larger than necessary, however. In a study of a 2-dimensional network [11], 20×20 links were large enough to document validity of the ergodic hypothesis. The more regular the system is, the smaller need the REV be. When the porous material consisted of a face centric array of large solid spheres, we found that half a unit cell was large enough to provide a REV [15]. The pores, or the substructure of the REV, might well be small. When the substructure consists of nanosized pores, the excess densities will depend on the curvature of the surfaces and of the three phase lines.

Additive contributions to the grand potential
Consider a REV of volume, V REV , with the temperature, T , the chemical potentials, µ n , µ w , µ r , kept constant by the surroundings. Additivity of the extensive variables can also be used to motivate the use of the grand potential. We use the grand canonical ensemble to calculate the grand canonical partition function, Ξ(T, V REV , µ n , µ w , µ r ) of the REV. The grand potential is given in terms of this partition function by The statistical mechanical treatment of small REV subsystems follows the procedure normal to large systems. The full distribution in phase space is normalized. The contributions to the REV grand potential are additive because the wetting and non-wetting fluids, the rock, the surfaces and the contact line are weakly coupled. Subsystems can be said to be weakly coupled when their interaction energy is small compared to the energies of the separate subsystems. The probability distribution in phase space is then equal to a product of the probability distributions in the subsystems For a given distribution of the fluids the coordinates of the subsystems are restricted to their volumes, surface areas or contact lines. It would then follow that the grand canonical partition function of the REV is the product of the partition functions of the subsystems. The origin of the configurational contributions is a probability distribution in phase space which also contains a factor that give the probability distribution of the fluids over all possible choices of subvolumes. By integrating over phase space, one also integrates over these distributions, which lead to an additional partition function. All of this together gives ×Ξ nr (T, Ω nr,REV , µ n , µ r )Ξ wr (T, Ω wr,REV , µ w , µ r )Ξ nw (T, Ω nw,REV , µ n , µ w ) As the masses of the components have no configurational contribution, it follows that Ξ conf does not depend on the chemical potentials.
When we now take the logarithm of Eq.24, multiply with −k B T , and identify the contributions of the subsystems with the corresponding contributions to the grand potential, we obtain Eq.16: +X nr,GC,REV (T, Ω nr,REV , µ n , µ r ) + X wr,GC,REV (T, Ω wr,REV , µ w , µ r ) +X nw,GC,REV (T, Ω nw,REV , µ n , µ w ) + X nwr,GC,REV (T, Λ nwr,REV , µ n , µ w , µ r ) We have indicated the variables, that the various contributions depend on. In the present case, we used that the components in the different bulk phases are immiscible. If they would be miscible, all material contributions depend on all the chemical potentials. Equations 17-20 are a direct consequence of Eq.25. The autonomous nature of the surfaces and three-phase contact lines is a direct consequence of the weak coupling between subsystems. Examples of weakly coupled systems can be found everywhere. The catalytic surface with its own temperature, the liquid vapor interface at steady state, are two examples [3].
The volume and the masses of the components only have material contributions. This is not a consequence of weak coupling, however. The analysis in the next section will show that the Gibbs energy has a material contribution only. All other thermodynamic energies, like the internal energy, the Helmholtz energy, and the enthalpy, have configurational contributions.

Gibbs equation. The meaning of equilibrium
The REV that we have considered so far is always in equilibrium with its surroundings. Away from global equilibrium the neighboring REVs are not necessarily in equilibrium with each other. The assumption of local equilibrium [3,4] covers this situation. It implies thermodynamic equilibrium within all REVs, and that the variables of the REV satisfy the Gibbs equation: dU GC,REV = T dS GC,REV − pdV REV + µ n dM n,GC,REV + µ w dM w,GC,REV + µ r dM r,GC,REV Each (additive) variable is constructed as described earlier. With this construction we can proceed as normal, for the sum of variables. We defend that Gibbs equation apply, using the argument that the equation applies for the single contributions and therefore also to the sum of additive variables. Given that the REV is large in Hill's sense U GC,REV is an Euler homogeneous function of the first order, leading to: By introducing the densities, Eq.26 becomes du GC = T ds GC + µ n dρ n,GC + µ w dρ w,GC + µ r dρ r,GC The corresponding Gibbs-Duhem relation is dp = s GC dT + ρ n,GC dµ n + ρ w,GC dµ w + ρ r,GC dµ r In the analysis we have argued that V It follows that G GC,REV has material contributions only. The coarse-grained variables used in this section, are those which will enter the thermodynamic description on the macro-level. They are relevant in equations of transport, like the Darcy equation, the Washburn equation, or equations that follow from a nonequilibrium thermodynamic expression, see below.
We need assume that the Gibbs equation is valid for the REV also when transport takes place. This is called the assumption of local equilibrium. Droplets can form at high flow rates, while ganglia may occur at low rates. There is a minimum size of the REV, for which Gibbs equation can be written [15]. When we assume that 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 [24]. Experimental evidence for the assumption was found by Erpelding [25].

Entropy production in porous media
We seen above how we can defined coarse-grained variables that describe a REV on the macroscale. Explicit expressions were given in terms of additive contributions. We did this for a REV in global equilibrium meaning that also T , V REV , µ n , µ w , µ r are well defined. The Gibbs equation 28 can then be written on the local form, in terms of densities. The discussion below considers n components. The time rate of change of the internal energy density is then from Gibbs equation Gradients in mass-and energy densities produce changes in the variables on the macro-scale. These lead to transport of heat and mass. The aim is to find the equations that govern this transport across the REV. The REV is chosen large enough to be macroscopic. The analysis follows a similar analysis for microporous media closely [14]. The balance equations for masses and internal energy densities of a REV are The Here J s is the entropy flux. Furthermore σ is the entropy productions. The entropy production is positive definite, σ ≥ 0 (the second law of thermodynamics).
We can now proceed to derive expressions for σ in the standard way [1,3], by combining the balance equations with Gibbs' equation from the previous section. We introduce the balance equations for mass and energy into Gibbs equation, see [1,3] for details. By comparing the result with the entropy balance, Eq.31, we identify first the entropy flux At the same time, we find the entropy production. Depending on the choice of independent variables, we can formulate the entropy production in various equivalent ways: The entropy production, σ, is here expressed using the energy flux, J u , as a variable, or using the measurable heat flux, J q as a variable. The final choice is motivated by practical wishes; what is measurable or computable. When we choose J u as flux, the conjugate force is ∂(1/T )/∂x, and the mass fluxes J i are driven by minus the gradient in the Planck potential, µ i /T . When, on the other hand we choose J q as the flux with the conjugate force ∂(1/T )/∂x, the mass fluxes are driven by minus the gradient in the chemical potential at constant temperature, divided by this temperature. The entropy production defines the independent thermodynamic driving forces and their conjugate fluxes. We have given two possible choices to demonstrate the flexibility [1,3,14]. The last expression is preferred for analysis of experiments.

Constitutive equations
From the entropy production we obtain the following flux-force relations when we use J u : According to Onsager [27], the conductivity matrix is symmetric: This is so, independent of the mechanism of transport. The magnitude of the coefficients reflect, of course, the mechanism of transport, but can not be used to conclude on the mechanism in play. When we use J q , we obtain This conductivity matrix is also symmetric: By using the relation between the J u and J q we can express one conductivity matrix, say eqs.36, in terms of the other matrix (of Eq. 38), or vice versa.

Fluctuation-dissipation theorems
We have so far have discussed how fluids flow through a porous medium, governed by the entropy production. All fluxes, densities, temperatures, pressures and chemical potentials were defined on the course-grained scale (not on the pore-scale). The description is now extended to incorporate fluctuations. We do this by adding fluctuating contributions to the fluxes, which are indicated by subscript R. These fluctuating contributions are crucial in the formulation of fluctuation-dissipation theorems [19][20][21]. Their are the level the quickly changing molecular contributions. Their temporal and spacial correlations are therefore on a molecular scale. Course grained this gives δ-functions in space and time on the macroscopic space and time scales. Their strength is given by 2k B times the Onsager conductivity matrix.
For the internal energy flux, we obtain The fluctuations and therefore the correlation of a flux pair, are characteristic for the pair. The correlation of fluctuation in the mass flux J r,R with the internal energy flux, differ from the corresponding correlation with the measurable heat flux. Self-correlations give information of diagonal coefficients, while crosscorrelations give information of the cross-coefficients. One set of fluctuation-dissipation relations can be derived from the other, and vice versa. The fluctuating flux-flux correlation functions decay fast compared to the macroscopic time scale [19][20][21]. The expressions apply for any condition, steady state or not, and do not refer to the nature of the flux-force relationships. The theorems apply, not only in equilibrium, but also away from equilibrium [1,4].

The chemical potential at constant temperature
The derivative of the chemical potential at constant temperature is needed in the driving forces in the previous section. For convenience we repeat its relation to the full chemical potential [1]. 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 entropies and volumes are equal to: and the last term of Eq. 46 is denoted by By reshuffling, we have the quantity of interest as the differential of the full chemical potential plus an entropic term; The Gibbs-Duhem's equation is 0 = SdT − V dp + n j=1 ρ j dµ j . By introducing Eq.49 into this equation we obtain an equivalent expression, to be used below:

The entropy production
Consider the case of two immiscible fluids, one more wetting (w) and one more non-wetting (n). The entropy production in Eq.35 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 the dimension of a velocity (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.50 to obtain 0 = ρ n dµ c n + ρ w dµ c w (53) The entropy production is invariant to the choice of variables. We can introduce the relations above and the explicit expression for dµ i,T into Eq.51, and find: where we used Eq.53 and the difference velocity J D : The difference 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.51 or 54, dictate the constitutive equations of the system.
The total volume flows through the REV Q n ≡ AJ n V n and Q w ≡ AJ w V w were used by Hansen et al [26]. Here A ≡ V / is the , and is the length of the REV in the flow direction. Hansen et al. [26] assumed that the total volume flow was an Euler homogeneous function of first order in the fractional areas, represent the saturation of the non-wetting and wetting fluid, respectively. Furthermore the porosity is given by φ = V p /V = A p /A. We will also use the total measurable heat flow Q q ≡ AJ q . 1 Because of the REV being homogeneous in the y, z directions, the total flows are the integrals of the flows across the cross-sectional area of the REV normal to the flow direction. By introducing the total flows in the expression for the total entropy production, Eq.51, we obtain: Similarly Eq.54 becomes where Q V ≡ AJ V and Q D ≡ AJ D . All the total fluxes as well as the corresponding thermodynamic forces only depend on x and t and not on y and z.

Constitutive equations
The entropy production given by Eqs.51 or 54, offer three equivalent choices for constitutive equations. For the sake of completeness, we give all of them.
Equation 51 dictates the following flux-force relations These are the equations given in Eq.38 for the special case of two fluids. According to Onsager [27,28], the conductivity matrix is symmetric. Equation 54 dictates the following flux-force relations Also this conductivity matrix is symmetric. Equation 56 dictates the following flux-force relations and Eq.57 gives Also these conductivity matrices are symmetric. From the relations between the fluxes, we find relations between the conductivity matrices of Eqs.58-61. In this way we find from Eqs.58 and 60 that and from Eqs.59 and 61 that The averages of the fluctuating contributions is zero. Their second moments satisfy the fluctuationdissipation theorem In Eqs.68 and 69 i, j are either q, w or n. The fluctuation-dissipation relations can be derived from Eq.65 by integrating over the cross section.
The averages of the fluctuating contributions is zero. Their second moments satisfy the fluctuationdissipation theorem In Eqs.70 and 71 i, j are either q, V or D. The fluctuation-dissipation relations can be derived from Eq.67 by integrating over the cross section.

Discussion and Conclusion
We have presented a definition of the representative elementary volume (REV) in terms of coarsegrained variables. For this set of variables, we have written the Gibbs equation for the REV, building on earlier work [14]. The additive nature of the contributions of the various phases, surface areas, contact lines and configurational contributions, was used to defend the procedure. The ergodic hypothesis applies. From Gibbs equation and the balance laws, we derived the entropy production for transport of heat and mass in porous media, following the standard procedure in non-equilibrium thermodynamics for heterogeneous systems [1][2][3].
On the basis provided by non-equilibrium thermodynamics, we formulated fluctuation dissipation theorems [14]. The theorems apply for a set of coarse-grained variables for the REV which obey Gibbs equation on the macro-scale. The theorems were specified for immiscible, non-isothermal, two-phase flow in a non-deformable medium.
These formulations of the fluctuation-dissipation theorems for porous media are new. We are thus not able to refer to supporting experimental work. It should be possible to determine the random fluctuations via optical techniques in a Hele-Shaw cell, however.
The very first support for these ideas was obtained from network simulations. Winkler et al. [29] simulated self-correlation and cross-correlation coefficients in a network with two-phase flow in a honeycomb lattice. The pore flow was modeled with the Washburn equation. The authors found that the matrix of coefficients obtained from auto and cross correlations of the fluctuating component flows was symmetric. This is the first result that indicates that Onsager reciprocal relations are obeyed for flow in a porous medium. The fluctuating flux-flux correlation functions in their study decayed fast on a macroscopic time scale (around 10 −3 s). In this paper we have shown that their integrals resulted in 2k B times the Onsager conductivities L P ij where i, j are either w or n. More systematic studies of this sort would be very instructive for the next steps to be taken. A crucial test will be to obtain the Onsager relations from the fluctuation-dissipation theorem, and from experiments defined directly from the flux-force relationships.
Other coarse-graining procedures exist. McClure et al. [23] and Khanamiri et al. [22] proposed to use a geometric state function for two-fluid flow in porous media based on Minkowski functionals, in order to characterize the complex arrangements of fluids and solid phases within a porous medium. The functionals were presented as separate, independent variables. Central variables are the volume, surface area and surface curvatures. To linear order, their variable construction seems similar to the one that follows from Hill's method [17,18]. In the traditional view, there is a viscous transport equation where pressure drop over connected phases drives the flow. In our approach, the REV pressure has contributions from several phases and interfaces, not only one.
Validity of Gibbs equation of the particular form used in this paper is essential, as a basis for all following derivations. In the different approach of McClure et al. [32], the Gibbs equation was formulated using time derivatives of the averages of the extensive measures, as well as fluctuation terms that arise due to the microscopic deviations in the intensive quantities. On this background, flux equations were written for the fluctuations themselves.
The L ij -coefficients used in this work are local coefficients. This means that they apply to a particular REV, and they are functions of the REV state variables that were defined earlier. According to Onsager, they are not functions of the driving forces. When these functions are known, it becomes possible to integrate the flux equation over a series of REVs. We may then arrive in a situation where the flux becomes a non-linear function of the overall driving force, as for chemical reactions. This property is not an argument against using relations which are linear on a local level, however.
The question has been raised in the literature about the origin of the energy dissipation as heat. In the present formulation, the dissipation is uniquely defined by the temperature of the surroundings, T 0 , times the entropy production [3]. The sources of the dissipation are then the product pairs of conjugate driving forces and fluxes. Each pair contributes to the dissipation. A typical example is the replacement of one fluid by another. The shift induces fluctuations, adding to the cross-correlation in Eq.65.
We have discussed that experiments and simulations can be done to test the proposed relations for consistency and performance. The advantage of constitutive equations presented here is that they can be used to determine permeabilities by experiments, and as well by the fluctuation-dissipation theorems for the flows. The two determinations should give the same result for the same overall driving force (pressure difference). Coefficients are often determined at steady state conditions. The equations do not depend on there being a steady state, however, as in the proposal of McClure et al. [32]. Once the coefficients are known, they can be used to model flow evolution in time on the time scale that is appropriate for the fluctuating flows.