A statistical mechanics for immiscible and incompressible two-phase flow in porous media

We construct a statistical mechanics for immiscible and incompressible two-phase flow in porous media under local steady-state conditions based on the Jaynes maximum entropy principle. A cluster entropy is assigned to our lack of knowledge of, and control over, the fluid and flow configurations in the pore space. As a consequence, two new variables describing the flow emerge: The agiture, that describes the level of agitation of the two fluids, and the flow derivative which is conjugate to the saturation. Agiture and flow derivative are the analogs of temperature and chemical potential in standard (thermal) statistical mechanics. The associated thermodynamics-like formalism reveals a number of hitherto unknown relations between the variables that describe the flow, including fluctuations. The formalism opens for new approaches to characterize porous media with respect to multi-phase flow for practical applications, replacing the simplistic relative permeability theory while still keeping the number of variables tractable.


I. INTRODUCTION
Flow in porous media [1][2][3][4] is a large field that spans many disciplines, from biology and chemistry to soil science, geophysics, materials science.A central problem is to find a theoretical description of immiscible two-phase flow in porous media at scales large enough for the pores to be negligible in size, often referred to as the continuum or Darcy scale.This is neither a new problem, nor are attempts at solutions new.The earliest attempt at solving the problem was that of Wyckoff and Botset [5] who regarded the flow of each of the immiscible fluids as one moving in a pore space reduced by the other fluids, thus reducing its own ability to move.This approach, now known as relative permeability theory, is today the standard framework used for practical applications.
There has been no lack of attempts to go beyond this simple theory.Perhaps the most advanced attempt to date is Thermodynamically Constrained Averaging Theory (TCAT) [6][7][8][9][10], based on thermodynamically consistent definitions made at the continuum scale based on volume averages of pore-scale thermodynamic quantities.Closure relations are then formulated at the macro-scale along the lines of the homogenization approach of Whit-taker [11].A key strength of TCAT is that all variables are defined in terms of pore-scale variables.A key disadvantage of TCAT is that many averaged variables are produced, and many complicated assumptions are needed to derive useful results.Another development based on non-equilibrium thermodynamics uses Euler homogeneity to define the up-scaled pressure.From this, Kjelstrup et al. derive constitutive equations for the flow while keeping the number of variables down [12][13][14].A challenge here is how to incorporate the structure of the fluid clusters spanning many pores.There is also a class of theories based on detailed and specific assumptions concerning the physics involved.An example is Local Porosity Theory [15][16][17][18][19][20].Another example is the Decomposition in Prototype Flow (DeProf) theory which is a fluid mechanical model combined with non-equilibrium statistical mechanics based on a classification scheme of fluid configurations at the pore level [21][22][23].
Recent work [24][25][26] explores a new approach to immiscible two-phase flow in porous media based on the Euler homogeneity theorem.It provides a transformation from the seepage velocity (i.e., average pore velocity) of the more wetting fluid, v w and the less wetting fluid v n to another pair of fluid velocities, the average seepage velocity of both fluids combined v p , and the co-moving velocity v m .The co-moving velocity appears as a result of the Euler homogeneity assumption.The transformation is reversible: knowing the average seepage velocity and the co-moving velocity, one can determine the seepage velocity of each fluid, (v p , v m ) (v w , v n ).
The mapping from the average velocity and the co-arXiv:2205.13791v1[physics.flu-dyn]27 May 2022 moving velocity to the seepage velocity of each fluid together with the assumption that the fluids are incompressible, constitutes a closed set of equations when supplemented by constitutive equations for the average velocity v p and the co-moving velocity v m .These two constitutive equations relate the two velocities to the driving forces, the pressure and saturation gradients.The constitutive equation for the average flow velocity v p relates to recent findings starting with Tallakstad et al. [27,28] who reported pressure drop proportional to the volumetric flow rate raised to a power 0.54 ± 0.08 in a two-dimensional model porous medium using a mixture of a compressible and an incompressible fluid under steady-state conditions.Aursjø et al. [29] found a power law with a somewhat larger exponent using the same model porous medium, but with two incompressible fluids.Similar results have since been observed by a number of groups, see [30][31][32].There has also been a considerable effort to understand these results theoretically and reproduce them numerically [27,28,[33][34][35][36][37][38][39][40][41].
The constitutive equation for the co-moving velocity v m has recently been studied by Roy et al. [26] by reverseengineering published relative permeability data and by using a dynamic pore network model.It turns out that this constitutive equation has a surprisingly simple form.We will return to this further on in this paper.
We note that these two constitutive laws are based on the collective behavior of both fluids combined.It would seem to be an impossible task to disassemble these two collective constitutive equations into one for each fluid.However, this is precisely what the mapping (v p , v m ) → (v w , v n ) makes possible.
The theory just described is modeled on thermodynamics as presented by Callen [42,43], which essentially consists of two ingredients: 1. an energy budget, i.e. the Gibbs relation and 2. an assumption that the variables describing the energy are Euler homogeneous functions.In the theory of Hansen et al. [24], only the second assumption is used.Hence, the following question may be posed: Can one complete the analogy between immiscible fluid flow in porous media and thermodynamics?Finding a positive answer to this question is the aim of this paper.
The theory we are developing here assumes local steady-state flow.This corresponds to local equilibrium in thermodynamics.By steady-state flow we mean that the macroscopic averages characterizing the flow are welldefined.This does not rule out that fluid clusters move, breaking up and merging at the pore level.
Finding an analogy between thermodynamics and nonthermal macroscopic systems is not new.Edwards and Oakeshott in their "Theory of Powders" [44] had the same aim when setting up a theory for static granular media.By making the assumption that all packing configurations having the same packing fraction are equally probable, thus constructing a microcanonical ensemble with the packing fraction playing the role of energy, an analogy with statistical mechanics was made.As a re-sult, a non-thermal pseudo-thermodynamics would follow.The ensuing theory thus had the same structure as thermodynamics, but the variables have nothing to do with thermodynamics -only the framework being the same.
In order to arrive at an analogue to thermodynamics in the present case, we take as a starting point the information theoretical statistical mechanics of Jaynes [45].This generalizes the principles of statistical mechanics from being a theory specifically for mechanical systems, e.g., assemblies of molecules, to a framework that can be implemented once certain criteria are in place, whatever the system.In essence the Jaynes approach generalizes the Laplace Principle of Insufficient Reason.Quoting Jaynes [45], "Laplace's 'Principle of Insufficient Reason' was an attempt to supply a criterion of choice, in which one said that two events are to be assigned equal probabilities if there is no reason to think otherwise." Jaynes furthermore builds his approach on Shannon's interpretation of entropy as a quantitative measure of what one does not know about a system -the less is known, the larger the entropy.From a set of properties that such a function of ignorance must have, Shannon managed to construct a unique one fulfilling them [46].One of these properties was that knowing nothing at all, this function would have its maximum when all possible states of the system would be equally probable, thus generalizing the Laplace principle of insufficient reason.Jaynes then interprets our knowledge about the system as constraints on our ignorance, so that that the Shannon entropy should be maximized under these constraints.
Constructing a theory along these principles for immiscible two-phase flow in porous media leads to a pseudothermodynamics describing the flow.This entails finding powerful relations between the variables describing the flow.Furthermore, it introduces new variables.One such new variable is the cluster entropy associated with the shapes of the fluid clusters.This allows us to define an agiture, essentially measuring the level of agitation of the two fluids.We have chosen the name "agiture" to emphasize that it is not a temperature.Another variable is the flow derivative which is the conjugate of the saturation.This variable is an analogue of the chemical potential in ordinary thermodynamics.A third new variable is the flow pressure which is conjugate to the porosity.
Furthermore, statistical mechanics is more powerful than thermodynamics in that all of thermodynamics may be derived from statistical mechanics but not vice versa.This is also the case here.For example, fluctuations in the macroscopic variables are accessible via statistical mechanics only.
We review in the next section the Euler homogeneity approach of [24].Here we describe the system we consider, defining central concepts.In particular, we derive the two-way mapping (v p , v m ) (v w , v n ) and discuss the co-moving velocity v m .In Section III we construct the Jaynes statistical mechanics for immiscible flow in porous media, and from this a pseudo-thermodynamics.In order to effectuate these ideas in practice, we need to define how we measure the relevant variables.This means defining averaging procedures, which should be done both in space and in time as already pointed out by McClure et al. [47,48].We go as far in this section as to define and exemplify the equivalents of the Maxwell relations in ordinary thermodynamics.The next section IV concerns saturation and porosity fluctuations.In Section V, we discuss the relation between the agiture, flow derivative and the pressure gradient.We do this by considering the internal balance in a porous medium having two regions with different matrix properties.
We end the paper by a discussion and conclusion section.Here we list a number of questions that remain open.

II. REVIEW OF EULER HOMOGENEITY APPROACH
Before we turn to constructing the statistical mechanics, we review the Euler homogeneity approach to immiscible and incompressible two-phase flow in porous media first introduced by Hansen et al. [24].

A. Representative Elementary Area
We imagine a porous medium plug as shown in Figure 1.It is homogeneous, i.e., the local porosity and permeability fluctuate around well-defined averages.The sides of the plug are sealed while the two ends remain open.We ignore gravity.A mixture of two immiscible fluids are injected through the lower end, and fluids are drained at the upper end.Flow into the plug is constant, leading to steady-state flow inside the plug [49].
Due to the sealed sides of the plug, the average flow direction is along the symmetry axis of the plug.We now imagine a plane cutting though the plug orthogonally to the average flow direction.In this plane, we choose a point.Around this point, we choose an area A, e.g., bounded by a circle, as shown in Figure 1.We assume that the area is large enough for averages of variables characterizing the flow are well-defined, but not larger.Furthermore, we assume the linear size of the area to be larger than any relevant correlation length in the system.This defines the Representative Elementary Area (REA) [50] at the chosen point.We may do the same at any point in the plane, and we may do this at any other such plane.

B. Some definitions
The REA has an area A. Part of this area is covered by the matrix material, whereas another part is covered by the pores.This latter area, A p , we will refer to as the FIG. 1: A Representative Elementary Area (REA) inside a porous medium plug.There is a volumetric flow rate Qp passing though the REA which has a transverse pore area Ap, of which an area Aw is covered by the wetting fluid.
transverse pore area.The porosity of the REA is then given by There is a volumetric flow rate through the REA shown in Figure 1, Q p .The seepage velocity of the fluids passing through the REA is then The flow consists of a mixture of two incompressible fluids, one being more wetting with respect to the matrix than the other one.We will refer to this fluid as the "wetting fluid."The less wetting fluid, we will refer to as the "non-wetting fluid."Each of them is associated with a volumetric flow rate, Q w and Q n , and we have The transverse pore area of the REA may also be divide into an area associated with the wetting fluid, A w and an area associated with the non-wetting fluid, A n , so that We may define a seepage velocity for each of the two fluids, We define the saturations of the two fluids passing through the REA as We may now combine these equations ( 2) -( 8) to expressing the average seepage velocity as We note that could have made these definitions based on the density of the two fluids, ρ w and ρ n rather than volumes.
The porosity φ, saturations S w and S n , and the seepage velocities v p , v w and v n are variables that may be associated with the chosen point surrounded by the REA.Since there for every point in the porous medium one may associate an REA, these variables, which do not depend on the size or shape of the REA, may be seen as continuous fields.
The variables we have defined should be well-defined.This is only possible if they are not fluctuating too strongly.Flooding processes typically generate fractal structures [4].There are, however, length scales associated with the mechanisms controlling these structures [51], and beyond the largest of these length scales, they cease to be fractal.The same is true for steady-state flow which we consider here.We expect self averaging to take place beyond these length scales, i.e., the relative strength of the fluctuations shrinks with increasing size of the system.Hence, we assume the plug and the REAs to be large enough for the fluctuations not to dominate.
In the discussion that follows, we need to assign another property to the variables beyond being self averaging.We need them to be state variables.By this we mean that they depend on the flow properties there and then and not the history of the flow.Erpelding et al. [49] studied this question experimentally and computationally, finding that this is indeed so.
A last aspect to be considered is that of hysteresis.There may indeed be hysteresis in the state variables.Knudsen and Hansen [52] studied immiscible two-phase flow under steady-state conditions using a dynamic pore network model.They found that there are two transitions between two-phase flow and single-phase flow when the wetting saturation is the control parameter.The transition between only the non-wetting fluid moving at low saturation to both fluids moving at higher saturation does not show any hysteresis with respect to which way one passes through the transition.On the other hand, the transition between only the wetting fluid moving at high saturation and both fluids moving at lower saturation does show a strong hysteresis, see Figure 2 in Knudsen and Hansen [52].This hysteresis is probably caused by this transition being equivalent to a first order or a spinodal phase transition.It is well known that such transitions in ordinary equilibrium thermodynamics lead to hysteretic behavior in the state variables.Hysteretic behavior signals that there are regions of parameter space where the macroscopic state variables are multi-valued.Hence, the underlying microscopic physics has more than one locally stable mode.

C. Relations between the seepage velocities
Hansen et al. [24] made the weak assumption that the volumetric flow rate through the REA is a homogeneous function of order one.That is, if we scale A → λA where λ is a scale factor, the volumetric flow rate would scale in the same way, i.e., We have here assumed that A w and A n are independent variables.This means that that we may change the area A of the REA by changing A w , while keeping A n fixed or changing A n while keeping A w fixed.This makes A and A p defined in equation ( 4) dependent variables.We refer to Hansen et al. [24] for details around this.By taking the derivative of equation ( 10) with respect to λ and then setting λ = 1, we find See Section 7.2 in Hansen et al. [24] for a step-by-step demonstration of how these derivatives are done for a capillary fiber bundle model.By invoking equations ( 2), ( 7) and ( 8), we find where we have defined We will refer to vw and vn as the thermodynamic velocities.
Let us now switch to treating S w and A p as the independent variables.Simple manipulations give Combining these two expressions with the definitions of the thermodynamic velocities ( 13) and ( 14) where we note that S n , vw , vn and v p are all function of S w and not of A p .
We combine equations ( 9) and ( 12), This does not imply that v w = vw and vn = v n .Rather, the most general relation between them is where v m is the co-moving velocity [24][25][26].
Combining these two expressions with equations ( 17) and ( 18) expresses the physical seepage velocities v w and v n in terms of the average seepage velocity v p and the co-moving velocity v m , We may invert equations ( 22) and ( 23), finding These four equations, ( 22) to (25), constitute a twoway mapping (v p , v m ) (v w , v n ).This means that having constitutive equations for v p and v m , makes it possible to derive constitutive equations for v w and v n .In other theories, such as relative permeability theory, only the mapping (v w , v n ) → v p is given, making it impossible to start from a constitutive equation for v p .
How to measure these variables from experimental data, see Roy et al. [26], who studied the co-moving velocity in detail using relative permeability data from the literature and from analyzing data obtained using a dynamic pore network model.They found that the comoving velocity takes the very simple form where A and B are coefficients dependent on the capillary number.It is still an open question as to why it has this simple form.We will in the following be able to give a partial answer.With these equations, the assumption that the two fluids are incompressible and constitutive equations for v p and v m , we have a closed set of equations describing the flow.

III. STATISTICAL MECHANICS A. Fluid configurations in a plug
We show in Figure 1 the porous plug that we discussed in Section II A. We will center our discussion on this plug in the following.
We assume that there is a volumetric flow rate Q p passing through the plug under steady-state conditions.This volumetric flow rate may be split into a wetting and a non-wetting volumetric flow rate, Q w and Q n so that Q p = Q w + Q n .We use script characters to distinguish the flow through the entire plug from the flow through an REA.We show in Figure 2 three planes orthogonal to the symmetry axis of the plug, i.e., orthogonal to the average flow direction.We introduce a coordinate z along the symmetry axis of the plug measuring the distance from the plug's lower boundary to any given plane.As the sides of the plug are sealed, we have that Q p is the same through the three orthogonal planes -or any other such orthogonal plane.The volumetric flow rate Q p is a conserved quantity as we move along the z axis.On the other hand, the volumetric flow rates of each of the two fluids, Q w and Q n are not conserved.Only their averages over several orthogonal planes will remain constant.One may imagine the fluids being layered in the flow direction to realize this.We also note that the transverse pore area A p will fluctuate from plane to plane due to fluctuations in the local porosity.The transverse area associated with the wetting fluid, A w , will also fluctuate for two reasons: A p fluctuates, but more importantly because the fluid clusters fluctuate.
We pick one of the planes, assuming it is at z 0 .We assume z 0 is large enough so that there are no end effects in the plane from the inlet at z = 0.One may at each point in the plane at z 0 assign a marker for whether the point is 1. in the solid matrix, 2. in the wetting fluid or 3. in the non-wetting fluid (thus ignoring the finite thickness of the interfaces and contact lines).We also assign a velocity vector to the point, which is zero if the point is in the solid matrix.This information collected for all points in the plane at z 0 at time t, we will refer to as the fluid configuration X = X (z 0 , t) in this plane.
We now imagine a stack of planes.We assume the neighboring planes are a distance dz > 0 apart.
In order to determine the flow configurations in each plane in the stack, it is not enough to know it at entry plane of the stack, z = z 0 .The reason for this is that we do not know the cluster structure inside the stack.However, the configurations in each plane are still measurable.
If we move through the plug along the z-axis where z 0 ≤ z ≤ z 1 , and the plug is long enough, we will explore the space of possible configurations X .We may do this at a fixed time t (hence moving infinitely fast) or at a finite speed along the z-axis so that time is running.We will in both cases explore the space of possible configurations.If we choose a given plane and then average over the configurations that pass it, we will not be averaging over the matrix structure, which will be fixed in the plane.However, if we imagine an ensemble of plugs each being a realization of the same statistical pore structure, we will eventually see all configurations X also in this case.
This leads us directly to defining a configurational probability density p(X ).We have also in the process sketched out an ergodicity assumption: the probability distribution for configurations in a given plane measured over an ensemble of plugs is the same as measuring it along different planes in a given plug.Time seems to have fallen out of the description.Time keeps track of the motion of each Lagrangian fluid element as we illustrate in Figure 2.This is more information than we need.All we need is to know the fluid configurations in the planes.For this purpose, the z-axis acts as an effective time axis.Hence, the reference to an equivalent to the ergodic hypothesis which normally is a statement about time averaging vs. ensemble averaging.
Since we cannot reconstruct the configurations inside the stack given a knowledge of the configuration at the entry plane at z 0 , one may be inclined to reject the idea to interpret the z-axis as an effective time direction.We point out that the same type of problem is encountered in relativistic mechanics of charged particles in strong fields where the position vs. time world lines of the particles are not single-valued functions of time [53,54].This means that it is not possible to reconstruct the later motion of such particles from knowing their configuration at a given time.
Statistical mechanics constitutes a calculational formalism based on the knowledge of the probability distribution for configurations.The Jaynes maximum entropy principle provides a recipe for determining the configurational probability density [45].Central to this principle is the definition of a function that quantitatively measures what we do not know about our system.This function is the Shannon entropy [46].We construct it for the present system as where the integral is over all hydrodynamically possible fluid configurations in the plane we focus on.The task is to calculate this entropy and from it determine p(X ).
We will refer to the entropy defined in equation ( 27) as the cluster entropy.We have chosen the name as it reflects the cluster structure that the fluids are making.We emphasize that the cluster entropy is not the thermodynamic entropy defined in other work such as [12][13][14].Whereas there is production of thermodynamic entropy in the plug as we are dealing with a driven system, there is no production of cluster entropy as we are dealing with steady-state flow.
The fluid states X are not discrete.Rather, they form a continuum.Hence, we use an integral in equation (27).There are mathematical problems related to defining the measure dX is this integral.However, these difficulties we presume are of the same type encountered in ordinary statistical mechanics.We will not delve into these problems here, but rather just note their existence and that they have been solved in ordinary statistical mechanics.

B. REA fluid configurations
We have in Section II A defined the REA.As for the entire plug, we may define a configuration X for the REA as a hydrodynamically possible configuration within the area A. The configuration in the plane where the REA sits is X .Hence, X is a subset of X .We also define the fluid configuration in the rest of the plane that is not part of the REA, X r .We will refer to this part of the plane as the "reservoir."Hence, we have that A central question is now, how independent are the configurations X and X r ?If they are independent, we may focus entirely on the REA configurations X as we may write the configurational probability for the entire plane as p(X ) = p(X)p r (X r ) , where p(X) is the configurational probability for the REA and p r (X r ) is the configurational probability for the reservoir.Fyhn et al. [55] have recently studied the validity of equation ( 29) in a two-dimensional dynamic pore network model.By changing the size of the two-dimensional plug, while keeping the size of the REA fixed, they checked whether the statistical distributions of Q p and A w were dependent on the size of the plug.Only a very weak dependency was found which decrease with size.It is therefore realistic to assume equation ( 29) to be valid for large enough plugs and REAs.
We proposed in Section III A to view the average flow direction through the plug, the z-axis as a playing the role of a time axis.Averaging over time thus corresponds to averaging over a stack of REAs as shown in Figure 3.When we in the next section refer to averaging, it is averaging over this stack we mean.
We may treat the interactions between the REA and the reservoir in different ways.We may remove the REA stack from the plug and treat it as a closed-off system.This amounts to treating the REA as a plug on its own.z FIG.3: Averaging over a stack of REAs in the flow (z) direction.This ensures that the averaging is also over fluctuations in the pore structure.
We may leave it in the original plug, allowing it to freely interact with the reservoir.We may allow the stack to interact fully with the reservoir, but in such a way that the amount of wetting fluid is kept constant.It is not possible to implement such constraints experimentally nor computationally, but theoretically it is.The same goes for the transverse pore area: we may keep it constant theoretically, but not experimentally or computationally.Nevertheless, we will see examples of such ensembles in the following, as they correspond to different control parameters.

C. Jaynes maximum entropy approach
Following the Jaynes recipe, we need to maximize the cluster entropy in the REA, which is Σ = − dXp(X) ln p(X) , (30) under the constraints of what we know about the system.The probability density p(X) is defined in equation (29).
There are three variables that we measure, the volumetric flow rate through the REA, Q p and the transverse pore area covered by the wetting fluid, A w and the transverse pore area A p .All three of them are averages over fluctuating quantities when performed as shown in Figure 3, i.e., we average over a stack of REAs.
The average of the volumetric flow rate is where the Q p (X) is associated with fluid configuration X.Likewise, the average wetting area is given by where A w (X) is the wetting area associated with configuration X.The third variable we consider is the average transverse pore area All three variables Q p , A w and A p are extensive in the area of the REA, A. The aim now is to determine p(X) based on the knowledge of these variable averages.Following Jaynes, we use the Lagrangian multiplier technique.We start by constructing the Lagrangian which we then will maximize, with the aim to determine p(X).We assume that the three Lagrange multipliers λ Q , λ w and λ A are intensive in the area of the REA, A. The Lagrange multiplier Λ, on the other hand, is extensive in the area A. It follows that the cluster entropy Σ is extensive in A.
Necessary conditions for maximizing the Lagrangian (34) are One may ask why not additional information is included here, such as the average wetting flow rate, Q w ?The answer to this lies in the Euler theory [24]: there are three independent variables in the problem for which we may fix their averages.We choose here Q p , A w and A p .Other averages will then follow from the formalism we are about to develop.
Equation (35) gives and the normalization condition, equation (36) gives where we have defined the partition function Z = Z(λ Q , λ w , λ A ).
We are now in a position to calculate the cluster entropy combining equations ( 30) and ( 40), Σ = − p(X) ln p(X) We note that These three equations may be solved to give λ Q , λ w and λ A as functions of the three variables we know, Q p , A w and A p .They also make equation ( 42) a triple Legendre transformation, changing the control variables λ Q → Q p and λ w → A w and λ A → A p .Hence, the control variables of the flow entropy are Σ = Σ(Q p , A w , A p ): We define a new variable It plays the role corresponding to that of a free energy in ordinary thermodynamics.Our next step is to invert equation ( 42) so that Q p becomes a function of Σ rather the other way round.That is, we transform Σ(Q p , A w , A p ) to Q p (Σ, A w , A p ).Hence, we may write equation ( 42) or (46) as where we have also used equation (47).We see that Hence, we note that the following equation, which forms part of the right hand side of equation ( 48), constitutes a Legendre transformation, Here Q F corresponds to another free energy in ordinary thermodynamics.
We rewrite equation (48) as The left hand side of this equation constitutes a Legendre transform, as we have Hence, we have now transformed equation ( 48) to D. Agiture, flow derivative and flow pressure Let us now define three new intensive variables built from the Lagrange multipliers λ Q , λ w and λ A , The first one, θ, by its resemblance to temperature in ordinary statistical mechanics, we will name the agiture.The unit of the agiture is the same as volumetric flow rate.However, it is an intensive variable.
The second one, π, we will name the flow pressure.This variable is the conjugate of the flow area A p , The unit of π is inverse velocity.Referring to equation (1), we see that A p is a measure of the porosity φ and π is therefore a velocity variable conjugate to the porosity.The third variable, equation (57) we name the flow derivative.We will explain the name in the next section.As π, its unit is that of a velocity.It is the conjugate of the transversal wetting fluid area A w , playing a role similar to a chemical potential in ordinary statistical mechanics.Through equation ( 7), we see that the flow derivative is the conjugate of the wetting saturation S w .

E. Connection with Euler homogeneity approach
We need to define one more variable, where Combining this definition with equations ( 50) and ( 54) gives We use the fact that both Q N (Σ, µ, A p ) and Q p (Σ, A w , A p ) are homogeneous functions of order one in the extensive variables Σ, A w and A p , We now set λ = 1/A p and combine these two expressions with equation (62), finding where we also have used equation (7) and we have defined the cluster entropy density We now divide equation (65) by A p , noting that i.e., v N and v p are velocities.We recognize v p = Q p /A p from equation (2) as the average seepage velocity of the two fluids.Equation (65) then takes on the form Let us now go back to equation (59) and use scaling relation (64) to find This expression is the reason why we name µ the flow derivative.We now compare these two equations ( 69) and (70) to equation ( 21), which we reproduce here: This equation is one of the central results derived by Hansen et al. [24].By comparison, we identify where the thermodynamic velocity of the non-wetting fluid, vn , is defined in equation (14).Equation ( 69) may be written and we have that These two equations are our central result.It demonstrates that the thermodynamic velocity of the nonwetting fluid is the Legendre transform of the average seepage velocity with respect to the wetting saturation and that the wetting saturation is minus the derivative of the non-wetting thermodynamic velocity with respect to the flow derivative.Equation ( 21) was derived based on the volumetric flow rate being a homogeneous function of the first kind in the transverse pore area.Now, we see this equation as a fundamental equation resulting from an underlying statistical mechanics, with variables (θ, Σ) and (µ, S w ), and (π, φ) forming conjugate pairs.
(74) We may split the integration over states X into an integral over A p and then over all states X that has a given where we have defined and where δ(• • • ) is the Dirac delta-function.We have used that Aφ = A p , see equation (1).We have that where Q M (θ, µ, A p ) is defined in equation ( 62).
We may now write partition function Z(θ, µ, A p ) in equation (76) as where we have defined which we may write as We may repeat this procedure one more time.We rewrite the partition function Z(θ, A w , A p ), equation (79) as where This microcanonical partition function (as Q p is the control variable) is also the (unnormalized) density of states with respect to the variables Q p , A w and A p .It demonstrates that all states X with the same Q p , A w and A p are equally probable.This brings us back to our starting point: the entropy has its maximum when all states that comply with the constraints (31), ( 32) and ( 33) are equally probable.

G. Co-moving velocity
The co-moving velocity, which is defined in equations ( 20) and (21), constitutes the bridge between the seepage velocities, equations ( 5) and ( 6), and the thermodynamics velocities, equations ( 13) and (14).By combining the defining equations for the co-moving velocities with the two equations ensuing from the Euler scaling assumption for the volumetric flow rate, we express the seepage velocity of the two fluid species in terms of the average seepage velocity and the co-moving velocity in equations ( 22) and (23).
Combining equation (23) with equation (73) based on statistical mechanics, we find a consistent structure when assuming that v n = v n (σ, µ) and v m = v m (σ, µ).Thus, we have expressed the physical seepage velocity for the non-wetting fluid within the pseudo-thermodynamic formalism we are developing.The corresponding physical seepage velocity for the wetting fluid may then be found e.g. by using equation (9).
The phenomenological form found by Roy et al. [26], equation (26), is consistent with the assumption for v m , as equation ( 26) then takes the form We note that the dependence of A and B on the capillary number is consistent with the two coefficients depending on the cluster entropy density σ.Why v m is linear in µ is not known.

H. Maxwell relations
We now see that the Euler homogeneity approach of Hansen et al. [24] was just the tip of an iceberg.Having anchored the approach in a statistical mechanics, we now have access to a rich formalism that parallels thermodynamics.
For example, we find the equivalents of the Maxwell relations in ordinary thermodynamics.We derive just one in the following.
We have that We combine this equation with equation (59), taking the cross derivatives, and finding which may be written as IV. FLUCTUATIONS AND AGITURE Our starting point are equations ( 41) and ( 47), i.e., = −θ ln dXe −Qp(X)/θ+µAw(X)/θ+πAp(X)/θ . ( This immediately gives Taking the derivative a second time with respect to µ gives where we have defined the fluctuations ∆A 2 w .We now combine equations ( 1) and ( 7) to find This allows us to transform equation (91) into We see that the agiture θ seems proportional to the fluctuations ∆(φS w ) 2 .However, due to the term (∂φS w /∂µ) θ,π , the relation between them is complex.We may calculate the porosity fluctuations by taking the partition function Z(θ, µ, π), equation (74), as starting point.We find where We note that the porosity field φ is a property of the matrix and not the flow.Hence, equation (94) gives a direct link between the agiture θ and the flow pressure π, θ = A∆φ 2 ∂π ∂φ θ,µ . (96)

V. CONDITIONS FOR STEADY STATE IN A HETEROGENEOUS PLUG
We consider in the following the conditions for steady state in a heterogeneous plug.
We show in Figure 4 a plug that is divided into a region A and a region B, which have different matrix properties.The difference may for example be that the wetting properties of the matrix in region A differ from those of region B. The full system, which consists of both regions A and B, is a closed system.
We will in the following focus on the four quantities Q p , A w , A p and Σ that describe the flow in the plug.Strictly speaking, we should change our notation, e.g., Q p → Q p , since we are considering the entire plug.However, in order to avoid complicating the notation and therefore making the material less accessible, we refrain from making the change.
These four quantities are extensive, i.e., additive.This means that we may express them in terms of the two regions A and B. We have We write Q p in terms of the control variables Σ, A w and A p , (101) Since the total volumetric flow rate is conserved in the plug, we must have that fluctuations in it must be zero, i.e, The fluctuations in Q A p and Q B p come from fluctuations in Σ A and Σ B , A A w and A B w , and A A p and A B p .We express δQ A p and δQ B p in terms of the fluctuations of these vari-ables, Likewise, we have There is no production of cluster entropy as this entropy is at a maximum.However, cluster entropy may move between regions A and B. Hence, we have δΣ = δΣ A + δΣ B = 0 .
(105) Furthermore, we keep the control variables A w and A p fixed.This is of course not possible experimentally.However, as a theoretical concept, it is permissible.Hence, we have that We now combine equations ( 102) -(107) to find We now assume that the cluster entropy fluctuations δΣ A are independent of the fluctuations in δA A w and δA A p .This leads to where we have used equations ( 49) and ( 55).
If we now furthermore assume that the fluctuations in δA A w and δA A p are uncorrelated, we find and These three criteria for steady state flow, equations (109), (110) and (111), are analogous to the criteria for two open systems in thermal contact and in equilibrium: here the temperature, pressure and the chemical potential must be the same.
It is important to point out the following.In ordinary thermodynamics, the rule is that the conjugate of a conserved quantity is constant in a heterogeneous system at equilibrium.This is e.g., the argument for the temperature being the same everywhere in the system at equilibrium as the entropy is at a maximum.This is the same argument as we have used here which leads to equation (109).However, consider two magnets with different magnetic susceptibilities in contact which is placed in a uniform magnetic field H.The conjugate of the magnetization M , which is the magnetic flux B, is not equal in the two magnets in contact, B A = B B .What goes wrong here is that it is not possible to let magnetization fluctuate between the two magnets so that its sum is constant, i.e., δM A = −δM B while simultaneously keeping the total internal energy and the total entropy constant, i.e., δU = δU A + δU B = 0 and δS = δS A + δS B = 0.In our system -steady-state two-phase flow in a porous medium -the following question then becomes central: is it possible to fulfill δQ p = 0 and δΣ = 0, and at the FIG. 4: The plug is now divided into a region A and a region B which differ in composition.The flow in region A is characterized by the variables θ A , µ A , π A and ∇P A , whereas in region B we have θ B , µ B , π B and ∇P B .Steady state flow is attained when θ A = θ B and µ A = µ B and π A = π B according to equations ( 109), ( 111) and ( 111).In addition we have same time have δA A w = −δA B w and δA A p = −δA B p ?The answer to this question is not obvious.Only experiments or computations on models may provide an answer.If the answer is no, only equation (109) will be valid, and not equations ( 110) and ( 111).
What about the pressure gradient?We must have that for steady-state conditions to apply.If the two pressure gradients were different, a pressure gradient orthogonal to the flow direction would develop, which would then divert the flow from the z direction, breaking the assumption of steady-state flow.Equation ( 112) is fulfilled if we have that due to (109), and if (110) and (111) are valid.

VI. DISCUSSION AND CONCLUSION
A central problem in the physics of porous media is how to scale up knowledge of the physics at the pore level to the continuum scale, often referred to as the Darcy scale, where the pores are negligible in size.The notion of up-scaling is of course only meaningful if we know which coarse-scale physics we are up-scaling to, and the traditional relative permeability theory has well known weak points.The theory of Hansen et al. [24] is a radically different approach to the coarse-scale physics.This theory is formally similar to some parts of thermodynamics, but it is the flow rates, and the relations between flow rates, that are the central players.In its core the theory is based on the volumetric flow rate being an Euler homogeneous function.
In the present paper, we have developed a statistical mechanics framework for immiscible incompressible two-phase flow.In this framework, total flow rate plays the same role as energy in ordinary statistical mechanics, and local steady state flow plays the role of local thermodynamic equilibrium.We have shown that the pseudo-thermodynamics that follows from the statistical mechanics framework extends the earlier theory of Hansen et al. [24], rendering it a complete pseudothermodynamics in the spirit of Edwards and Oakeshot's pseudo-thermodynamics for powders [44].
In the same way that ordinary statistical mechanics serves as a tool for calculating macro-scale properties from known molecular scale physics, the present statistical mechanics links the pore scale hydrodynamic description to a continuum scale physics, thus solving the up-scaling problem.The key link is the partition function Z(θ, µ, π) defined in equation (41).The integral is over all physical fluid configurations, and he physics sits in the measure dX over configurations and pore structure.
The up-scaling from the intermittent flow of fluid clusters through pore space to a description with a small number of continuum scale variables admits an immense level of ignorance.This ignorance is reflected in the pseudo-thermodynamics as the cluster entropy.We call the conjugate variable to this cluster entropy the agiture.The agiture measure the level of agitation in the flow.Typically a high agiture is associated with high flow rates, in the same way high temperature is associated with high energies in thermal systems.
We have not proposed here any technique as to how the agiture θ may be measured, or controlled, experimentally or computationally.The flow pressure π seems also difficult to measure.Measuring the flow derivative is, on the other hand, seemingly easier to measure.In terms of the average seepage velocity and the saturation, it may be written µ = ∂v p (σ, S w ) ∂S w σ . (114) Since we are assuming steady-state flow, the cluster entropy is at a maximum and therefore constant.We also note that we have used five different ensembles in this work, where control parameters are respectively (θ, µ, π) , (θ, µ, A p ) , (Σ, µ, A p ) , (θ, A w , A p ) , (Σ, A w , A p ) .
The first of these ensembles is easily realizable experimentally as this describes an REA communicating freely with the rest of the porous medium.The second one, with A p controlled, is feasible e.g. by using 3D printing techniques to construct the porous medium.The two last three ensembles must be seen as theoretical constructs.We note, however, that this is standard procedure in thermodynamics in that one always starts with the Gibbs relation, which relates variations in the internal energy to variations in the extensive variables irrespective of whether this ensemble is accessible experimentally or not.
We see that a steady-state condition in a heterogeneous porous medium such as that shown in Figure 4 is that the agiture is constant throughout the entire system, equation (109).On the other hand, the two additional equilibrium conditions, equations (110) and (111), hinge on additional assumptions that need to be verified experimentally or computationally.In light of the above discussion, verifying the validity of equation ( 110) is the easier one of the three equations.
The probability to find a flow configuration p(X) in the REA which is the central to the theory we present here.It was emphasized at the beginning of Section III C that the fluid configuration X refers to the configuration in the REA, not the entire plane cutting through the plug, X = X ∪ X r .It is crucial for the theory that p(X) does not depend on the properties of the entire plug.Recent numerical work by Fyhn et al. [55] indicates that this is indeed correct.
Given all these caveats and open question, the fact remains that the framework we have presented here, opens up for viewing immiscible two-phase flow in porous media in a different way than before.We have here divided the problem into 1.constructing a framework by which the concepts that are necessary are put in place, and 2. connecting these concepts to the underlying fluid dynamics and thermodynamics.This paper accomplishes the first goal.

FIG. 2 :
FIG.2:We show three planes orthogonal to the symmetry axis of the porous plug from Figure1at different z coordinates.We also show the trajectories of two Lagrangian fluid elements.They started simultaneously at time t = 0 at the plug inlet.At time t they are at the indicated positions.