Advances in Water Resources 171 (2023) 104336

Contents lists available at ScienceDirect

Advances in Water Resources

journal homepage: www.elsevier.com/locate/advwatres

A statistical mechanics framework for immiscible and incompressible

two-phase flow in porous media
Alex Hansen a ,∗, Eirik Grude Flekkøy b , Santanu Sinha b , Per Arne Slotte c
a
PoreLab, Department of Physics, NTNU, N-7491 Trondheim, Norway
b
PoreLab, Department of Physics, University of Oslo, N-0316 Oslo, Norway
c
PoreLab, Department of Geoscience and Petroleum, NTNU, N-7491 Trondheim, Norway

ARTICLE INFO ABSTRACT

Keywords: We construct a statistical mechanics for immiscible and incompressible two-phase flow in porous media under
Multiphase flow in porous media local steady-state conditions based on the Jaynes maximum entropy principle. A cluster entropy is assigned
Jaynes statistical mechanics to our lack of knowledge of, and control over, the fluid and flow configurations in the pore space. As a
Cluster entropy
consequence, two new variables describing the flow emerge: The agiture, which describes the level of agitation
Agiture
of the two fluids, and the flow derivative, which is conjugate to the saturation. Agiture and flow derivative are
Thermodynamic description
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.

1. Introduction variables. A key disadvantage of TCAT is that many averaged variables

are produced, and many complicated assumptions are needed to derive
Flow in porous media (Bear, 1988; Sahimi, 2011; Blunt, 2017; useful results.
Feder et al., 2022) is a large field that spans many disciplines, from Barenblatt et al. (2002) point out that the key assumption in relative
biology and chemistry to soil science, geophysics, materials science. permeability theory is that the flow is locally in a steady state, even
A central problem is to find a theoretical description of immiscible if the flow as a whole is developing. This allows the central variables
two-phase flow in porous media at scales large enough for the pores
of that theory, the relative permeability and the capillary pressure, to
to be negligible in size, often referred to as the continuum or Darcy
be functions of the saturation alone. They then go on to generalize the
scale. This is neither a new problem, nor are attempts at solutions
theory to flow which is locally in out of equilibrium, exploring how the
new. The earliest attempt at solving the problem was that of Wyckoff
and Botset (1936) who regarded the flow of each of the immiscible central variables change. Wang et al. (2019) take these ideas further by
fluids as one moving in a pore space reduced by the other fluids, thus introducing dynamic length scales due to the mixing zone variations,
reducing its own ability to move. This approach, now known as relative over which the spatial averaging is done.
permeability theory, is today the standard framework used for practical Another development based on non-equilibrium thermodynamics
applications. uses Euler homogeneity to define the up-scaled pressure. From this,
There has been no lack of attempts to go beyond this simple theory. Kjelstrup et al. derive constitutive equations for the flow while keeping
One of the most advanced attempt to date is Thermodynamically the number of variables down (Kjelstrup et al., 2018, 2019; Bedeaux
Constrained Averaging Theory (TCAT) (Hassanizadeh and Gray, 1990, and Kjelstrup, 2022). A challenge here is how to incorporate the
1993a,b; Niessner et al., 2011; Gray and Miller, 2014), based on structure of the fluid clusters spanning many pores.
thermodynamically consistent definitions made at the continuum scale There is also an ongoing effort in constructing a scaled-up theory
based on volume averages of pore-scale thermodynamic quantities.
based on geometrical properties, or more precisely the Minkowski
Closure relations are then formulated at the macro-scale along the
functionals (McClure et al., 2018; Armstrong et al., 2019; McClure
lines of the homogenization approach of Whitaker (1986). A key
et al.).
strength of TCAT is that all variables are defined in terms of pore-scale

∗ Corresponding author.
E-mail addresses: [email protected] (A. Hansen), [email protected] (E.G. Flekkøy), [email protected] (S. Sinha), [email protected] (P.A. Slotte).

https://doi.org/10.1016/j.advwatres.2022.104336
Received 20 June 2022; Received in revised form 11 October 2022; Accepted 26 October 2022
Available online 2 November 2022
0309-1708/© 2022 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

There is also a class of theories based on detailed and specific aim when setting up a theory for static granular media. By making the
assumptions concerning the physics involved. An example is Local assumption that all packing configurations having the same packing
Porosity Theory (Hilfer and Besserer, 2000; Hilfer, 2006,a,b; Hilfer fraction are equally probable, thus constructing a microcanonical en-
and Döster, 2010; Döster et al., 2012). Another example is the De- semble with the packing fraction playing the role of energy, an analogy
composition in Prototype Flow (DeProf) theory which is a fluid me- with statistical mechanics was made. As a result, a non-thermal pseudo-
chanical model combined with non-equilibrium statistical mechanics thermodynamics would follow. The ensuing theory thus had the same
based on a classification scheme of fluid configurations at the pore structure as thermodynamics, but the variables have nothing to do with
level (Valavanides et al., 1998; Valavanides, 2012, 2018). thermodynamics — only the framework being the same.
Recent work (Hansen et al., 2018; Roy et al., 2020, 2022) explores In order to arrive at an analogue to thermodynamics in the present
a new approach to immiscible two-phase flow in porous media based case, we take as a starting point the information theoretical statistical
on the Euler homogeneity theorem. It provides a transformation from mechanics of Jaynes (1957). This generalizes the principles of statistical
the seepage velocity (i.e., average pore velocity) of the more wetting mechanics from being a theory specifically for mechanical systems,
fluid, 𝑣𝑤 and the less wetting fluid 𝑣𝑛 to another pair of fluid velocities, e.g., assemblies of molecules, to a framework that can be implemented
the average seepage velocity of both fluids combined 𝑣𝑝 , and the co- once certain criteria are in place, whatever the system. In essence
moving velocity 𝑣𝑚 . The co-moving velocity appears as a result of the Jaynes approach generalizes the Laplace Principle of Insufficient
the Euler homogeneity assumption. The transformation is reversible: Reason. Quoting Jaynes (1957), ‘‘Laplace’s ‘Principle of Insufficient
knowing the average seepage velocity and the co-moving velocity, one Reason’ was an attempt to supply a criterion of choice, in which one
can determine the seepage velocity of each fluid, (𝑣𝑝 , 𝑣𝑚 ) ⇄ (𝑣𝑤 , 𝑣𝑛 ). said that two events are to be assigned equal probabilities if there is no
The mapping from the average velocity and the co-moving velocity reason to think otherwise’’.
to the seepage velocity of each fluid together with the assumption that Jaynes furthermore builds his approach on Shannon’s interpretation
the fluids are incompressible, constitutes a closed set of equations when of entropy as a quantitative measure of what one does not know about a
supplemented by constitutive equations for the average velocity 𝑣𝑝 system — the less is known, the larger the entropy. From a set of prop-
and the co-moving velocity 𝑣𝑚 . These two constitutive equations relate erties that such a function of ignorance must have, Shannon managed to
the two velocities to the driving forces, the pressure and saturation construct a unique one fulfilling them (Shannon, 1948). One of these
gradients. properties was that knowing nothing at all, this function would have
The constitutive equation for the average flow velocity 𝑣𝑝 relates to its maximum when all possible states of the system would be equally
recent findings starting with Tallakstad et al. (2009a,b) who reported probable, thus generalizing the Laplace principle of insufficient reason.
pressure drop proportional to the volumetric flow rate raised to a power Jaynes then interprets our knowledge about the system as constraints
0.54 ± 0.08 in a two-dimensional model porous medium using a mix- on our ignorance, so that the Shannon entropy should be maximized
ture of a compressible and an incompressible fluid under steady-state under these constraints.
conditions. Aursjø et al. (2014) found a power law with a somewhat Constructing a theory along these principles for immiscible two-
larger exponent using the same model porous medium, but with two phase flow in porous media leads to a pseudo-thermodynamics de-
incompressible fluids. Similar results have since been observed by a scribing the flow. This entails finding powerful relations between the
number of groups, see Sinha et al. (2017), Gao et al. (2020), Zhang variables describing the flow. Furthermore, it introduces new variables.
et al. (2021). There has also been a considerable effort to understand One such new variable is the cluster entropy associated with the shapes
these results theoretically and reproduce them numerically (Tallakstad of the fluid clusters. This allows us to define an agiture, essentially
et al., 2009a,b; Grøva and Hansen, 2011; Sinha and Hansen, 2012; measuring the level of agitation of the two fluids. We have chosen the
Sinha et al., 2013; Xu and Wang, 2014; Yiotis et al., 2019; Roy et al., name ‘‘agiture’’ which is a contraction of the two words ‘‘agitation’’ and
2019; Fyhn et al., 2021; Lanza et al., 2022; Roy et al.; Cheon et al.). ‘‘temperature’’ to emphasize that it is not a temperature in the usual
The constitutive equation for the co-moving velocity 𝑣𝑚 has recently sense. Another variable is the flow derivative which is the conjugate of
been studied by Roy et al. (2022) by reverse-engineering published the saturation. This variable is an analogue of the chemical potential
relative permeability data and by using a dynamic pore network model. in ordinary thermodynamics. A third new variable is the flow pressure
It turns out that this constitutive equation has a surprisingly simple which is conjugate to the porosity.
form. We will return to this further on in this paper. Furthermore, statistical mechanics is more powerful than thermody-
We note that these two constitutive laws are based on the collective namics in that all of thermodynamics may be derived from statistical
behavior of both fluids combined. It would seem to be an impossible mechanics but not vice versa. This is also the case here. For example,
task to disassemble these two collective constitutive equations into one fluctuations in the macroscopic variables are accessible via statistical
for each fluid. However, this is precisely what the mapping (𝑣𝑝 , 𝑣𝑚 ) → mechanics only.
(𝑣𝑤 , 𝑣𝑛 ) makes possible. We review in the next section the Euler homogeneity approach
The theory just described is modeled on thermodynamics as pre- of Hansen et al. (2018). Here we describe the system we consider, defin-
sented by Callen (1974, 1991), which essentially consists of two ing central concepts. In particular, we derive the two-way mapping
ingredients: 1. an energy budget, i.e. the Gibbs relation and 2. an (𝑣𝑝 , 𝑣𝑚 ) ⇄ (𝑣𝑤 , 𝑣𝑛 ) and discuss the co-moving velocity 𝑣𝑚 . In Section 3
assumption that the variables describing the energy are Euler homoge- we construct the Jaynes statistical mechanics for immiscible flow in
neous functions. In the theory of Hansen et al. (2018), only the second porous media, and from this a pseudo-thermodynamics. In order to
assumption is used. Hence, the following question may be posed: Can effectuate these ideas in practice, we need to define how we mea-
one complete the analogy between immiscible fluid flow in porous sure the relevant variables. This means defining averaging procedures,
media and thermodynamics? Finding a positive answer to this question which should be done both in space and in time as already pointed out
is the aim of this paper. by McClure et al. (2021a,b). We go as far in this section as to define
The theory we are developing here assumes local steady-state flow. and exemplify the equivalents of the Maxwell relations in ordinary
This corresponds to local equilibrium in thermodynamics. By steady- thermodynamics. The next Section 4 concerns saturation and porosity
state flow we mean that the macroscopic averages characterizing the fluctuations. In Section 5, we discuss the relation between the agiture,
flow are well-defined. This does not rule out that fluid clusters move, flow derivative and the pressure gradient. We do this by considering the
break up and merge at the pore level. internal balance in a porous medium having two regions with different
Finding an analogy between thermodynamics and non-thermal matrix properties.
macroscopic systems is not new. Edwards and Oakeshott in their We end the paper by a discussion and conclusion section. Here we
‘‘Theory of Powders’’ (Edwards and Oakeshott, 1989) had the same list a number of questions that remain open.

2
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

The transverse pore area of the REA may also be divide into an area
associated with the wetting fluid, 𝐴𝑤 and an area associated with the
non-wetting fluid, 𝐴𝑛 , so that

𝐴𝑝 = 𝐴𝑤 + 𝐴𝑛 . (4)
We may define a seepage velocity for each of the two fluids,
𝑄𝑤
𝑣𝑤 = , (5)
𝐴𝑤
𝑄
𝑣𝑛 = 𝑛 . (6)
𝐴𝑛
Fig. 1. A Representative Elementary Area (REA) inside a porous medium plug. There We define the saturations of the two fluids passing through the REA as
is a volumetric flow rate 𝑄𝑝 passing though the REA which has a transverse pore area
𝐴𝑝 , of which an area 𝐴𝑤 is covered by the wetting fluid.
𝐴𝑤
𝑆𝑤 = , (7)
𝐴𝑝
𝐴
𝑆𝑛 = 𝑛 . (8)
2. Review of Euler homogeneity approach 𝐴𝑝
We may now combine these Eqs. (2) —(8) to express the average
Before we turn to constructing the statistical mechanics, we review seepage velocity as
the Euler homogeneity approach to immiscible and incompressible 𝑣𝑝 = 𝑆𝑤 𝑣𝑤 + 𝑆𝑛 𝑣𝑛 . (9)
two-phase flow in porous media first introduced by Hansen et al.
We note that could have made these definitions based on the density
(2018).
of the two fluids, 𝜌𝑤 and 𝜌𝑛 rather than volumes.
The porosity 𝜙, saturations 𝑆𝑤 and 𝑆𝑛 , and the seepage velocities
2.1. Representative elementary area 𝑣𝑝 , 𝑣𝑤 and 𝑣𝑛 are variables that may be associated with the chosen
point surrounded by the REA. Since there for every point in the porous
We imagine a porous medium plug as shown in Fig. 1. It is ho- medium one may associate an REA, these variables, which do not
mogeneous, i.e., the local porosity and permeability fluctuate around depend on the size or shape of the REA, may be seen as continuous
fields.
well-defined averages. The sides of the plug are sealed while the two
The variables we have defined should be well-defined. This is only
ends remain open. We ignore gravity. A mixture of two immiscible
possible if they not fluctuating too strongly, so that it is possible
fluids are injected through the lower end, and fluids are drained at the
to measure their averages. By this we mean that the average does
upper end. Flow into the plug is constant, leading to steady-state flow not depend on the number of samplings. Flooding processes typically
inside the plug (Erpelding et al., 2013). generate fractal structures (Feder et al., 2022). There are, however,
Due to the sealed sides of the plug, the average flow direction is length scales associated with the mechanisms controlling these struc-
along the symmetry axis of the plug. We now imagine a plane cutting tures (Måløy et al., 2021), and beyond the largest of these length scales,
though the plug orthogonally to the average flow direction. In this they cease to be fractal. The same is true for steady-state flow which
plane, we choose a point. Around this point, we choose an area 𝐴, we consider here. We expect self averaging to take place beyond these
e.g., bounded by a circle, as shown in Fig. 1. We assume that the area length scales, i.e., the relative strength of the fluctuations with respect
is large enough for averages of variables characterizing the flow are to their averages shrinks with increasing size of the system. Hence, we
well-defined, but not larger. Furthermore, we assume the linear size assume the plug and the REAs to be large enough for the fluctuations
of the area to be larger than any relevant correlation length in the not to dominate. That such a combination of sizes is possible to find, has
recently been investigated by Fyhn et al. using a dynamic pore network
system. This defines the Representative Elementary Area (REA) (Bear and
model.
Bachmat, 2012) at the chosen point. We may do the same at any point
In the discussion that follows, we need to assign another property
in the plane, and we may do this at any other such plane. 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
2.2. Some definitions and then and not the history of the flow. Erpelding et al. (2013) studied
this question experimentally and computationally, finding that this is
indeed so.
The REA has an area 𝐴. Part of this area is covered by the matrix
A last aspect to be considered is that of hysteresis. There may
material, whereas another part is covered by the pores. This latter area,
indeed be hysteresis in the state variables. Knudsen and Hansen (2006)
𝐴𝑝 , we will refer to as the transverse pore area. The porosity of the REA
studied immiscible two-phase flow under steady-state conditions us-
is then given by ing a dynamic pore network model. They found that there are two
𝐴𝑝 transitions between two-phase flow and single-phase flow when the
𝜙= . (1) wetting saturation is the control parameter. The transition between
𝐴
only the non-wetting fluid moving at low saturation to both fluids
There is a volumetric flow rate through the REA shown in Fig. 1,
moving at higher saturation does not show any hysteresis with respect
𝑄𝑝 . The seepage velocity of the fluids passing through the REA is then to which way one passes through the transition. On the other hand,
𝑄𝑝 the transition between only the wetting fluid moving at high satura-
𝑣𝑝 = . (2) tion and both fluids moving at lower saturation does show a strong
𝐴𝑝
hysteresis, see Figure 2 in Knudsen and Hansen (2006). This hysteresis
The flow consists of a mixture of two incompressible fluids, one is probably caused by this transition being equivalent to a first order
being more wetting with respect to the matrix than the other one. We or a spinodal phase transition. It is well known that such transitions
will refer to this fluid as the ‘‘wetting fluid’’. The less wetting fluid, we in ordinary equilibrium thermodynamics lead to hysteretic behavior in
will refer to as the ‘‘non-wetting fluid’’. Each of them is associated with the state variables. Hysteretic behavior signals that there are regions
a volumetric flow rate, 𝑄𝑤 and 𝑄𝑛 , and we have of parameter space where the macroscopic state variables are multi-
valued. Hence, the underlying microscopic physics has more than one
𝑄𝑝 = 𝑄𝑤 + 𝑄𝑛 . (3) locally stable mode.

3
A. Hansen et al. Advances in Water Resources 171 (2023) 104336
( )
𝑑𝑣𝑝
2.3. Relations between the seepage velocities 𝑣𝑛 = 𝑣𝑝 − 𝑆 𝑤 − 𝑣𝑚 . (23)
𝑑𝑆𝑤
Hansen et al. (2018) made the weak assumption that the volumetric We may invert Eqs. (22) and (23), finding
flow rate through the REA is a homogeneous function of order one. That
is, if we scale 𝐴 → 𝜆𝐴 where 𝜆 is a scale factor, the volumetric flow 𝑣𝑝 = 𝑆𝑤 𝑣𝑤 + 𝑆𝑛 𝑣𝑛 , (24)
rate would scale in the same way, i.e., 𝑑𝑣 𝑑𝑣
𝑣𝑚 = 𝑆𝑤 𝑤 + 𝑆𝑛 𝑛 . (25)
𝑑𝑆𝑤 𝑑𝑆𝑤
𝑄𝑝 (𝜆𝐴𝑤 , 𝜆𝐴𝑛 ) = 𝜆𝑄𝑝 (𝐴𝑤 , 𝐴𝑛 ) . (10)
These four equations, (22) to (25), constitute a two-way mapping
We have here assumed that 𝐴𝑤 and 𝐴𝑛 are independent variables. (𝑣𝑝 , 𝑣𝑚 ) ⇄ (𝑣𝑤 , 𝑣𝑛 ). This means that having constitutive equations for
This means that we may change the area 𝐴 of the REA by changing 𝑣𝑝 and 𝑣𝑚 , makes it possible to derive constitutive equations for 𝑣𝑤
𝐴𝑤 , while keeping 𝐴𝑛 fixed or changing 𝐴𝑛 while keeping 𝐴𝑤 fixed. and 𝑣𝑛 . In other theories, such as relative permeability theory, only the
This makes 𝐴 and 𝐴𝑝 defined in Eq. (4) dependent variables. We refer mapping (𝑣𝑤 , 𝑣𝑛 ) → 𝑣𝑝 is given, making it impossible to start from a
to Hansen et al. (2018) for details around this. constitutive equation for 𝑣𝑝 .
By taking the derivative of Eq. (10) with respect to 𝜆 and then How to measure these variables from experimental data, see Roy
setting 𝜆 = 1, we find et al. (2022), who studied the co-moving velocity in detail using
( ) ( )
𝜕𝑄𝑝 𝜕𝑄𝑝 relative permeability data from the literature and from analyzing data
𝑄𝑝 (𝐴𝑤 , 𝐴𝑛 ) = 𝐴𝑤 + 𝐴 . (11)
𝜕𝐴𝑤 𝐴 𝜕𝐴𝑛 𝐴 𝑛 obtained using a dynamic pore network model. They found that the
𝑛 𝑤
co-moving velocity takes the very simple form
See Section 7.2 in Hansen et al. (2018) for a step-by-step demonstration
of how these derivatives are done for a capillary fiber bundle model. 𝑑𝑣𝑝
𝑣𝑚 = 𝐴 + 𝐵 , (26)
By invoking Eqs. (2), (7) and (8), we find 𝑑𝑆𝑤
( ) ( ) where 𝐴 and 𝐵 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.
= 𝑆𝑤 𝑣̂ 𝑤 + 𝑆𝑛 𝑣̂ 𝑛 , (12) With these equations, the assumption that the two fluids are incom-
where we have defined pressible and constitutive equations for 𝑣𝑝 and 𝑣𝑚 , we have a closed set
( ) of equations describing the flow.
𝜕𝑄𝑝
𝑣̂ 𝑤 = , (13)
𝜕𝐴𝑤 𝐴
( )
𝑛 3. Statistical mechanics
𝜕𝑄𝑝
𝑣̂ 𝑛 = . (14)
𝜕𝐴𝑛 𝐴 3.1. Fluid configurations in a plug
𝑤

We will refer to 𝑣̂ 𝑤 and 𝑣̂ 𝑛 as the thermodynamic velocities.

We now change variables from (𝐴𝑤 , 𝐴𝑛 ) being independent to We show in Fig. 1 the porous plug that we discussed in Section 2.1.
(𝑆𝑤 , 𝐴𝑝 ) being the independent variables. We then have We will center our discussion on this plug in the following.
𝐴𝑤 (𝑆𝑤 , 𝐴𝑝 ) = 𝑆𝑤 𝐴𝑝 and 𝐴𝑛 (𝑆𝑤 , 𝐴𝑝 ) = (1 − 𝑆𝑊 )𝐴𝑝 . We find We assume that there is a volumetric flow rate 𝑝 passing through
( ) ( ) ( ) the plug under steady-state conditions. This volumetric flow rate may
𝜕 𝜕 𝑆 𝜕
= + 𝑛 , (15) be split into a wetting and a non-wetting volumetric flow rate, 𝑤
𝜕𝐴𝑤 𝐴 𝜕𝐴𝑝 𝑆 𝐴𝑝 𝜕𝑆𝑤 𝐴
( )
𝑛
( ) 𝑤 ( )
𝑝 and 𝑛 so that 𝑝 = 𝑤 + 𝑛 . We use script characters to distinguish
𝜕 𝜕 𝑆 𝜕 the flow through the entire plug from the flow through an REA. We
= − 𝑤 . (16)
𝜕𝐴𝑛 𝐴 𝜕𝐴𝑝 𝑆 𝐴𝑝 𝜕𝑆𝑤 𝐴 show in Fig. 2 three planes orthogonal to the symmetry axis of the
𝑤 𝑤 𝑝
plug, i.e., orthogonal to the average flow direction. We introduce
Combining these two expressions with the definitions of the thermody-
a coordinate 𝑧 along the symmetry axis of the plug measuring the
namic velocities (13) and (14) gives
distance from the plug’s lower boundary to any given plane. As the
𝑑𝑣𝑝 sides of the plug are sealed, we have that 𝑝 is the same through the
𝑣̂ 𝑤 = 𝑣𝑝 + 𝑆𝑛 , (17)
𝑑𝑆𝑤 three orthogonal planes — or any other such orthogonal plane. The
𝑑𝑣𝑝 volumetric flow rate 𝑝 is a conserved quantity as we move along the
𝑣̂ 𝑛 = 𝑣𝑝 − 𝑆𝑤 , (18)
𝑑𝑆𝑤 𝑧 axis. On the other hand, the volumetric flow rates of each of the two
where we note that 𝑆𝑛 , 𝑣̂ 𝑤 , 𝑣̂ 𝑛 and 𝑣𝑝 are all function of 𝑆𝑤 and not of fluids, 𝑤 and 𝑛 are not conserved. Only their averages over several
𝐴𝑝 . orthogonal planes will remain constant. One may imagine the fluids
We combine Eqs. (9) and (12), being layered in the flow direction to realize this. We also note that
the transverse pore area 𝐴𝑝 will fluctuate from plane to plane due to
𝑣𝑝 = 𝑆𝑤 𝑣𝑤 + 𝑆𝑛 𝑣𝑛 = 𝑆𝑤 𝑣̂ 𝑤 + 𝑆𝑛 𝑣̂ 𝑛 . (19) fluctuations in the local porosity. The transverse area associated with
the wetting fluid, 𝐴𝑤 , will also fluctuate for two reasons: 𝐴𝑝 fluctuates,
This does not imply that 𝑣𝑤 = 𝑣̂ 𝑤 and 𝑣̂ 𝑛 = 𝑣𝑛 . Rather, the most general
but more importantly because the fluid clusters fluctuate.
relation between them is
We pick one of the planes at 𝑧 = 𝑧0 . We introduce a two-dimensional
𝑣̂ 𝑤 = 𝑣𝑤 + 𝑣𝑚 𝑆𝑛 , (20) grid that divides the plane into voxels. Each voxel is assigned a set of
𝑣̂ 𝑛 = 𝑣𝑛 − 𝑣𝑚 𝑆𝑤 , (21) parameters, its area, its transverse pore area, the area covered by the
wetting fluid and volumetric flow rate through it. The values assigned
where 𝑣𝑚 is the co-moving velocity (Hansen et al., 2018; Roy et al., 2020, to the voxels, which are coarse grained at the scale of the voxel, are
2022). the result of physical fluid configurations: the respect the underlying
Combining these two expressions with Eqs. (17) and (18) expresses hydrodynamics and thermodynamics. For example, these values may be
the physical seepage velocities 𝑣𝑤 and 𝑣𝑛 in terms of the average seepage the values assigned to the nodes in a Lattice Boltzmann model system
velocity 𝑣𝑝 and the co-moving velocity 𝑣𝑚 , or it may be the status of the links in a dynamic pore network model.
( ) It may also be the status of the pixels in transversal CT scan through a
𝑑𝑣𝑝
𝑣𝑤 = 𝑣𝑝 + 𝑆 𝑛 − 𝑣𝑚 , (22) porous plug.
𝑑𝑆𝑤

4
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

This information assigned to each voxels in the plane at 𝑧0 at time

𝑡, we will refer to as the fluid configuration  = (𝑧0 , 𝑡).
We now imagine a stack of planes. We assume the neighboring
planes are a distance 𝑑𝑧 > 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, 𝑧 = 𝑧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 𝑧-axis where 𝑧0 ≤ 𝑧 ≤ 𝑧1 , and
the plug is long enough, we will explore the space of possible config-
urations . We may do this at a fixed time 𝑡 (hence moving infinitely
fast) or at a finite speed along the 𝑧-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 Fig. 2. We show three planes orthogonal to the symmetry axis of the porous plug from
the plane. However, if we imagine an ensemble of plugs each being a Fig. 1 at different 𝑧 coordinates. We also show the trajectories of two Lagrangian fluid
realization of the same statistical pore structure, we will eventually see elements. They started simultaneously at time 𝑡 = 0 at the plug inlet. At time 𝑡 they
are at the indicated positions.
all configurations  also in this case.
This leads us directly to defining a configurational probability den-
sity 𝑝().
̃ 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 Fig. 2.
This is more information than we need. All we need is to know the
fluid configurations in the planes. For this purpose, the 𝑧-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 𝑧0 ,
one may be inclined to reject the idea to interpret the 𝑧-axis as an Fig. 3. Averaging over a stack of REAs in the flow (𝑧) direction. This ensures that the
effective time direction. We point out that the same type of problem averaging is also over fluctuations in the pore structure.
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 (Feynman, 1948; Hansen and Ravndal, 3.2. REA fluid configurations
1981). This means that it is not possible to reconstruct the later motion
of such particles from knowing their configuration at a given time. We have in Section 2.1 defined the REA. As for the entire plug,
Statistical mechanics constitutes a calculational formalism based on we may define a configuration 𝑋 for the REA as a hydrodynamically
the knowledge of the probability distribution for configurations. The possible configuration within the area 𝐴. The configuration in the plane
Jaynes maximum entropy principle provides a recipe for determining where the REA sits is . Hence, 𝑋 is a subset of . We also define the
the configurational probability density (Jaynes, 1957). Central to this fluid configuration in the rest of the plane that is not part of the REA,
principle is the definition of a function that quantitatively measures
𝑋𝑟 . We will refer to this part of the plane as the ‘‘reservoir’’. Hence, we
what we do not know about our system. This function is the Shannon
have that
entropy (Shannon, 1948). We construct it for the present system as
 = 𝑋 ∪ 𝑋𝑟 . (28)
𝛴𝑝𝑙𝑎𝑛𝑒 = − 𝑑 𝑝()
̃ ln 𝑝()
̃ , (27)
∫
A central question is now, how independent are the configurations 𝑋
where the integral is over all hydrodynamically possible fluid configu- and 𝑋𝑟 ? If they are independent, we may focus entirely on the REA
rations in the plane we focus on. The task is to calculate this entropy configurations 𝑋 as we may write the configurational probability for
and from it determine 𝑝().
̃ the entire plane as
We will refer to the entropy defined in Eq. (27) as the cluster entropy.
We have chosen the name as it reflects the cluster structure that the 𝑝()
̃ = 𝑝(𝑋)𝑝𝑟 (𝑋𝑟 ) , (29)
fluids are making. We emphasize that the cluster entropy is not the
where 𝑝(𝑋) is the configurational probability for the REA and 𝑝𝑟 (𝑋𝑟 ) is
thermodynamic entropy defined in other work such as (Kjelstrup et al.,
the configurational probability for the reservoir.
2018, 2019; Bedeaux and Kjelstrup, 2022). Whereas there is production
of thermodynamic entropy in the plug as we are dealing with a driven Fyhn et al. have recently studied the validity of Eq. (29) in a two-
system, there is no production of cluster entropy as we are dealing with dimensional dynamic pore network model. By changing the size of
steady-state flow. the two-dimensional plug, while keeping the size of the REA fixed,
The fluid states  are not discrete. Rather, they form a continuum. they checked whether the statistical distributions of 𝑄𝑝 and 𝐴𝑤 were
Hence, we use an integral in Eq. (27). There are mathematical problems dependent on the size of the plug. Only a very weak dependency
related to defining the measure 𝑑 is this integral. However, these was found which decrease with size. It is therefore realistic to assume
difficulties we presume are of the same type encountered in ordinary Eq. (29) to be valid for large enough plugs and REAs.
statistical mechanics. We will not delve into these problems here, but We proposed in Section 3.1 to view the average flow direction
rather just note their existence and that they have been solved in through the plug, the 𝑧-axis as a playing the role of a time axis.
ordinary statistical mechanics. Averaging over time thus corresponds to averaging over a stack of REAs

5
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

𝜕
as shown in Fig. 3. When we in the next section refer to averaging, it = 0, (36)
𝜕𝛬
is averaging over this stack we mean. 𝜕
We may treat the interactions between the REA and the reservoir = 0, (37)
𝜕𝜆𝑄
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 = 0, (38)
𝜕𝜆𝑤
plug on its own. We may leave it in the original plug, allowing it to 𝜕
freely interact with the reservoir. We may allow the stack to interact = 0. (39)
𝜕𝜆𝐴
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 One may ask why not additional information is included here,
experimentally nor computationally, but theoretically it is. The same such as the average wetting flow rate, 𝑄𝑤 ? The answer to this lies in
goes for the transverse pore area: we may keep it constant theoretically, the Euler theory (Hansen et al., 2018): there are three independent
but not experimentally or computationally. Nevertheless, we will see variables in the problem for which we may fix their averages. We
examples of such ensembles in the following, as they correspond to choose here 𝑄𝑝 , 𝐴𝑤 and 𝐴𝑝 . Other averages will then follow from the
different control parameters. formalism we are about to develop. The thermodynamic formalism we
therefore are about to develop will then concern the thermodynamic
3.3. Jaynes maximum entropy approach velocities, Eqs. (13) and (14). We will then need the co-moving velocity
𝑣𝑚 to translate the results into the physical velocities.
Following the Jaynes recipe, we need to maximize the cluster Eq. (35) gives
entropy in the REA, which is 𝑝(𝑋) = 𝑒1−𝛬 𝑒−𝜆𝑄 𝑄𝑝 (𝑋)−𝜆𝑤 𝐴𝑤 (𝑋)−𝜆𝐴 𝐴𝑝 (𝑋) , (40)

𝛴=− 𝑑𝑋𝑝(𝑋) ln 𝑝(𝑋) , (30) and the normalization condition, Eq. (36) gives
∫
under the constraints of what we know about the system. The proba- 𝑒𝛬−1 = 𝑑𝑋𝑒−𝜆𝑄 𝑄𝑝 (𝑋)−𝜆𝑤 𝐴𝑤 (𝑋)−𝜆𝐴 𝐴𝑝 (𝑋) = 𝑍 , (41)
∫
bility density 𝑝(𝑋) is defined in Eq. (29).
There are three variables that we measure, the volumetric flow where we have defined the partition function 𝑍 = 𝑍(𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ).
rate through the REA, 𝑄𝑝 and the transverse pore area covered by the We are now in a position to calculate the cluster entropy combining
wetting fluid, 𝐴𝑤 and the transverse pore area 𝐴𝑝 . All three of them Eqs. (30) and (40),
are averages over fluctuating quantities when performed as shown in
Fig. 3, i.e., we average over a stack of REAs. 𝛴 = − 𝑝(𝑋) ln 𝑝(𝑋)
∫
The average of the volumetric flow rate is
= ln 𝑍(𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) + 𝜆𝑄 𝑄𝑝 + 𝜆𝑤 𝐴𝑤 + 𝜆𝐴 𝐴𝑝 , (42)
𝑑𝑋 𝑝(𝑋)𝑄𝑝 (𝑋) = 𝑄𝑝 . (31) We note that
∫
( )
where the 𝑄𝑝 (𝑋) is associated with fluid configuration 𝑋. Likewise, the 𝜕 ln 𝑍
𝑄𝑝 (𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) = − , (43)
average wetting area is given by 𝜕𝜆𝑄 𝜆 ,𝜆
( ) 𝑤 𝐴

𝜕 ln 𝑍
𝑑𝑋 𝑝(𝑋)𝐴𝑤 (𝑋) = 𝐴𝑤 , (32) 𝐴𝑤 (𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) = − , (44)
∫ 𝜕𝜆𝑤 𝜆 ,𝜆
𝑄 𝐴
( )
where 𝐴𝑤 (𝑋) is the wetting area associated with configuration 𝑋. The 𝜕 ln 𝑍
𝐴𝑝 (𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) = − . (45)
third variable we consider is the average transverse pore area 𝜕𝜆𝐴 𝜆 ,𝜆
𝑄 𝑤

𝑑𝑋 𝑝(𝑋)𝐴𝑝 (𝑋) = 𝐴𝑝 . (33) These three equations may be solved to give 𝜆𝑄 , 𝜆𝑤 and 𝜆𝐴 as functions
∫ of the three variables we know, 𝑄𝑝 , 𝐴𝑤 and 𝐴𝑝 . They also make Eq. (42)
All three variables 𝑄𝑝 , 𝐴𝑤 and 𝐴𝑝 are extensive in the area of the REA, a triple Legendre transformation, changing the control variables 𝜆𝑄 →
𝐴. The aim now is to determine 𝑝(𝑋) based on the knowledge of these 𝑄𝑝 and 𝜆𝑤 → 𝐴𝑤 and 𝜆𝐴 → 𝐴𝑝 . Hence, the control variables of the flow
variable averages. entropy are 𝛴 = 𝛴(𝑄𝑝 , 𝐴𝑤 , 𝐴𝑝 ):
Following Jaynes, we use the Lagrangian multiplier technique. We
start by constructing the Lagrangian which we then will maximize, 𝛴(𝑄𝑝 , 𝐴𝑤 , 𝐴𝑝 ) = − 𝑝(𝑋) ln 𝑝(𝑋)
∫
( )
= − 𝑑𝑋 𝑝(𝑋) ln 𝑝(𝑋) 𝜕 ln 𝑍
∫ = ln 𝑍(𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) − 𝜆𝑄
( ) 𝜕𝜆𝑄 𝜆 ,𝜆
( ) ( )𝑤 𝐴
− 𝛬 1− 𝑑𝑋 𝑝(𝑋) 𝜕 ln 𝑍 𝜕 ln 𝑍
∫ − 𝜆𝑤 − 𝜆𝐴 , (46)
( ) 𝜕𝜆𝑤 𝜆 ,𝜆
𝑄 𝐴
𝜕𝜆𝐴 𝜆 ,𝜆
𝑄 𝑤
+ 𝜆𝑄 𝑄𝑝 − 𝑑𝑋 𝑝(𝑋)𝑄𝑝 (𝑋)
∫ We define a new variable 𝑄𝐺 = 𝑄𝐺 (𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) through the equation
( )
+ 𝜆𝑤 𝐴 𝑤 − 𝑑𝑋 𝑝(𝑋)𝐴𝑤 (𝑋)
∫
( ) 𝑍(𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) = 𝑒−𝜆𝑄 𝑄𝐺 (𝜆𝑄 ,𝜆𝑤 ,𝜆𝐴 ) . (47)
+ 𝜆𝐴 𝐴 𝑝 − 𝑑𝑋 𝑝(𝑋)𝐴𝑝 (𝑋) , (34)
∫ It plays the role corresponding to that of a free energy in ordinary
thermodynamics.
with the aim to determine 𝑝(𝑋).
Our next step is to invert Eq. (42) so that 𝑄𝑝 becomes a function of
We assume that the three Lagrange multipliers 𝜆𝑄 , 𝜆𝑤 and 𝜆𝐴 are
𝛴 rather the other way round. That is, we transform 𝛴(𝑄𝑝 , 𝐴𝑤 , 𝐴𝑝 ) to
intensive in the area of the REA, 𝐴. The Lagrange multiplier 𝛬, on the
𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ). Hence, we may write Eq. (42) or (46) as
other hand, is extensive in the area 𝐴. It follows that the cluster entropy
𝛴 is extensive in 𝐴. 1
𝑄𝐺 (𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) = 𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ) − 𝛴
Necessary conditions for maximizing the Lagrangian (34) are 𝜆𝑄
𝜕 𝜆𝑤 𝜆
= 0, (35) + 𝐴𝑤 + 𝐴𝑝 𝐴 , (48)
𝜕𝑝(𝑋) 𝜆𝑄 𝜆𝑄

6
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

where we have also used Eq. (47). playing a role similar to a chemical potential in ordinary statistical
We see that mechanics. Through Eq. (7), we see that the flow derivative is the
( ) conjugate of the wetting saturation 𝑆𝑤 .
𝜕𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ) 1
= . (49)
𝜕𝛴 𝐴 ,𝐴 𝜆𝑄
𝑤 𝑝
3.5. Connection with Euler homogeneity approach
Hence, we note that the following equation, which forms part of the
right hand side of Eq. (48), constitutes a Legendre transformation, We need to define one more variable,
( ) 1 ( )
𝑄𝐹 𝜆𝑄 , 𝐴𝑤 , 𝐴𝑝 = 𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ) − 𝛴 𝑄𝑁 (𝛴, 𝜇, 𝐴𝑝 ) = 𝑄𝑀 𝜃, 𝜇, 𝐴𝑝 − 𝛴(𝜃, 𝜇, 𝐴𝑝 )𝜃 , (60)
𝜆𝑄
( )
𝜕𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ) where
= 𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ) − 𝛴 . (50) ( )
𝜕𝛴 𝐴 ,𝐴
𝜕𝑄𝑀 (𝜃, 𝜇, 𝐴𝑝 )
𝑤 𝑝 𝛴(𝜃, 𝜇, 𝐴𝑝 ) = . (61)
𝜕𝜃 𝜇,𝐴𝑝
Here 𝑄𝐹 corresponds to another free energy in ordinary thermodynam-
ics. Combining this definition with Eqs. (50) and (54) gives
We rewrite Eq. (48) as ( ) ( )
𝑄𝑁 𝛴, 𝜇, 𝐴𝑝 = 𝑄𝑝 𝛴, 𝐴𝑤 , 𝐴𝑝 − 𝐴𝑤 𝜇 . (62)
𝜆𝐴
𝑄𝐺 (𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) − 𝐴𝑝
𝜆𝑄 We use the fact that both 𝑄𝑁 (𝛴, 𝜇, 𝐴𝑝 ) and 𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ) are ho-
( ) 𝜆 mogeneous functions of order one in the extensive variables 𝛴, 𝐴𝑤 and
= 𝑄𝐹 𝜆𝑄 , 𝐴𝑤 , 𝐴𝑃 + 𝐴𝑤 𝑤 . (51)
𝜆𝑄 𝐴𝑝 ,
The left hand side of this equation constitutes a Legendre transform, ( ) ( )
𝑄𝑁 𝜆𝛴, 𝜇, 𝜆𝐴𝑝 = 𝜆𝑄𝑁 𝛴, 𝜇, 𝐴𝑝 , (63)
( ) ( )
( ) 𝐴𝑝 𝑄𝑝 𝜆𝛴, 𝜆𝐴𝑤 , 𝜆𝐴𝑝 = 𝜆𝑄𝑝 𝛴, 𝐴𝑤 , 𝐴𝑝 . (64)
𝑄𝑀 𝜆𝑄 , 𝜆𝑤 , 𝐴𝑝 = 𝑄𝐺 (𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) − 𝜆 , (52)
𝜆𝑄 𝐴 We now set 𝜆 = 1∕𝐴𝑝 and combine these two expressions with Eq. (62),
as we have finding
( )
⎛ ⎞ 𝐴𝑝 𝑄𝑁 (𝜎, 𝜇, 1) = 𝐴𝑝 𝑄𝑝 𝜎, 𝑆𝑤 , 1 − 𝐴𝑤 𝜇 , (65)
⎜ 𝜕𝑄𝐺 (𝜆𝑄 , 𝜆𝑤 , 𝜆𝐴 ) ⎟
⎜ ( ) ⎟ = 𝐴𝑝 . (53) where we also have used Eq. (7) and we have defined the cluster
𝜆
⎜ 𝜕 𝜆𝐴 ⎟
⎝ 𝑄 ⎠𝜆𝑄 ,𝜆𝑤 entropy density

Hence, we have now transformed Eq. (48) to 𝛴

𝜎= . (66)
𝐴𝑝
( ) ( ) 𝜆
𝑄𝑀 𝜆𝑄 , 𝜆𝑤 , 𝐴𝑝 = 𝑄𝐹 𝜆𝑄 , 𝐴𝑤 , 𝐴𝑝 + 𝐴𝑤 𝑤 . (54) We now divide Eq. (65) by 𝐴𝑝 , noting that
𝜆𝑄

3.4. Agiture, flow derivative and flow pressure 𝑄𝑁 (𝜎, 𝜇, 1) = 𝑣𝑁 (𝜎, 𝜇) , (67)
𝑄𝑝 (𝜎, 𝑆𝑤 , 1) = 𝑣𝑝 (𝜎, 𝑆𝑤 ) , (68)
Let us now define three new intensive variables built from the
Lagrange multipliers 𝜆𝑄 , 𝜆𝑤 and 𝜆𝐴 , i.e., 𝑣𝑁 and 𝑣𝑝 are velocities. We recognize 𝑣𝑝 = 𝑄𝑝 ∕𝐴𝑝 from Eq. (2) as
the average seepage velocity of the two fluids. Eq. (65) then takes on
1
𝜃 = + , (55) the form
𝜆𝑄
𝜆 𝑣𝑁 (𝜎, 𝜇) = 𝑣𝑝 (𝜎, 𝑆𝑤 ) − 𝑆𝑤 𝜇 . (69)
𝜋 = − 𝐴 , (56)
𝜆𝑄
Let us now go back to Eq. (59) and use scaling relation (64) to find
𝜆 ( )
𝜇 = − 𝑤 . (57) 𝜕𝑣𝑝 (𝜎, 𝑆𝑤 )
𝜆𝑄 =𝜇. (70)
𝜕𝑆𝑤 𝜎
The first one, 𝜃, by its resemblance to temperature in ordinary
statistical mechanics, we will name the agiture. The unit of the agiture This expression is the reason why we name 𝜇 the flow derivative. We
is the same as volumetric flow rate. However, it is an intensive variable. now compare these two Eqs. (69) and (70) to Eq. (21), which we
The second one, 𝜋, we will name the flow pressure. This variable is reproduce here:
( )
the conjugate of the flow area 𝐴𝑝 , 𝑑𝑣𝑝
( ) 𝑣̂ 𝑛 = 𝑣𝑝 − 𝑆𝑤 . (71)
𝜕𝑄𝐹 (𝜃, 𝐴𝑤 , 𝐴𝑝 ) 𝑑𝑆𝑤
𝜋 =
𝜕𝐴𝑝 𝜃,𝐴𝑤
This equation is one of the central results derived by Hansen et al.
( ) (2018). By comparison, we identify
𝜕𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 )
= . (58)
𝜕𝐴𝑝 𝛴,𝐴 𝑣𝑁 = 𝑣̂ 𝑛 , (72)
𝑤

The unit of 𝜋 is inverse velocity. Referring to Eq. (1), we see that 𝐴𝑝 where the thermodynamic velocity of the non-wetting fluid, 𝑣̂ 𝑛 , is
is a measure of the porosity 𝜙 and 𝜋 is therefore a velocity variable defined in Eq. (14). Eq. (69) may be written
conjugate to the porosity.
The third variable, Eq. (57) we name the flow derivative. We will 𝑣̂ 𝑛 (𝜎, 𝜇) = 𝑣𝑝 (𝜎, 𝑆𝑤 (𝜎, 𝜇)) − 𝑆𝑤 (𝜎, 𝜇)𝜇 , (73)
explain the name in the next section. As 𝜋, its unit is that of a velocity. and we have that
It is the conjugate of the transversal wetting fluid area 𝐴𝑤 , ( )
( ) 𝜕 𝑣̂ 𝑛 (𝜎, 𝜇)
𝜕𝑄𝐹 (𝜃, 𝐴𝑤 , 𝐴𝑝 ) 𝑆𝑤 (𝜎, 𝜇) = − . (74)
𝜇 = 𝜕𝜇 𝜎
𝜕𝐴𝑤 𝜃,𝐴𝑝 These two equations are our central result. It demonstrates that the ther-
( )
𝜕𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ) modynamic velocity of the non-wetting fluid is the Legendre transform of the
= . (59)
𝜕𝐴𝑤 𝛴,𝐴 average seepage velocity with respect to the wetting saturation and that the
𝑝

7
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

wetting saturation is minus the derivative of the non-wetting thermodynamic 3.7. Co-moving velocity
velocity with respect to the flow derivative. Eq. (21) was derived based on
the volumetric flow rate being a homogeneous function of the first kind The co-moving velocity, which is defined in Eqs. (20) and (21), con-
in the transverse pore area. Now, we see this equation as a fundamen- stitutes the bridge between the seepage velocities, Eqs. (5) and (6), and
tal equation resulting from an underlying statistical mechanics, with the thermodynamics velocities, Eqs. (13) and (14). By combining the
variables (𝜃, 𝛴) and (𝜇, 𝑆𝑤 ), and (𝜋, 𝜙) forming conjugate pairs. 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
3.6. Partition functions average seepage velocity and the co-moving velocity in Eqs. (22) and
(23).
Combining Eq. (23) with Eq. (74) based on statistical mechanics, we
The partition function is given by find a consistent structure
( ) ( )
𝑍(𝜃, 𝜇, 𝜋) = 𝑑𝑋 𝑒−𝑄𝑝 (𝑋)∕𝜃+𝜇𝐴𝑤 (𝑋)∕𝜃+𝜋𝐴𝑝 (𝑋)∕𝜃 . (75) 𝑣𝑛 (𝜎, 𝜇) = 𝑣𝑝 𝜎, 𝑆𝑤 (𝜎, 𝜇) − 𝑆𝑤 (𝜎, 𝜇) 𝜇 − 𝑣𝑚 (𝜎, 𝜇) , (84)
∫
We may split the integration over states 𝑋 into an integral over 𝐴𝑝 and when assuming that 𝑣𝑛 = 𝑣𝑛 (𝜎, 𝜇) and 𝑣𝑚 = 𝑣𝑚 (𝜎, 𝜇). Thus, we
have expressed the physical seepage velocity for the non-wetting fluid
then over all states 𝑋 that has a given 𝐴𝑝 ,
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 Eq. (9).
= 𝑑𝐴𝑝 𝑒𝜋𝐴𝑝 ∕𝜃 𝑍(𝜃, 𝜇, 𝐴𝑝 ) , (76) The phenomenological form found by Roy et al. (2022), Eq. (26), is
∫0
consistent with the assumption for 𝑣𝑚 , as Eq. (26) then takes the form
where we have defined
𝑣𝑚 (𝜎, 𝜇) = 𝐴(𝜎) + 𝐵(𝜎)𝜇 , (85)
𝑍(𝜃, 𝜇, 𝐴𝑝 )
with
( )
= 𝑑𝑋 𝛿(𝐴𝑝 (𝑋) − 𝐴𝑝 ) 𝑒−𝑄𝑝 (𝑋)∕𝜃+𝜇𝐴𝑤 (𝑋)∕𝜃 𝜕𝑣𝑚 (𝜎, 𝜇)
∫ 𝐵(𝜎) = . (86)
𝜕𝜇 𝜎
1 −𝑄𝑝 (𝑋)∕𝜃+𝜇𝐴𝑤 (𝑋)∕𝜃
= 𝑑𝑋 𝛿(𝜙(𝑋) − 𝜙) 𝑒 , (77)
𝐴∫ We note that the dependence of 𝐴 and 𝐵 on the capillary number is
consistent with the two coefficients depending on the cluster entropy
and where 𝛿(⋯) is the Dirac delta-function. We have used that 𝐴𝜙 = 𝐴𝑝 ,
density 𝜎.
see Eq. (1). We have that
Why 𝑣𝑚 is linear in 𝜇 is not known.
1 −𝑄𝑀 (𝜃,𝜇,𝐴𝑝 )∕𝜃
𝑍(𝜃, 𝜇, 𝐴𝑝 ) = 𝑒 , (78)
𝐴 3.8. Maxwell relations
where 𝑄𝑀 (𝜃, 𝜇, 𝐴𝑝 ) is defined in Eq. (62).
We may now write partition function 𝑍(𝜃, 𝜇, 𝐴𝑝 ) in Eq. (77) as We now see that the Euler homogeneity approach of Hansen et al.
(2018) was just the tip of an iceberg. Having anchored the approach
𝐴𝑝
in a statistical mechanics, we now have access to a rich formalism that
𝑍(𝜃, 𝜇, 𝐴𝑝 ) = 𝑑𝐴𝑤 𝑒𝜇𝐴𝑤 ∕𝜃 𝑍(𝜃, 𝐴𝑤 , 𝐴𝑝 ) , (79)
∫0 parallels thermodynamics.
For example, we find the equivalents of the Maxwell relations in
where we have defined
ordinary thermodynamics. We derive just one in the following.
𝑍(𝜃, 𝐴𝑤 , 𝐴𝑝 ) We have that
( )
𝜕𝑄𝐹 (𝜃, 𝐴𝑤 , 𝐴𝑝 )
= 𝑑𝑋 𝛿(𝐴𝑤 (𝑋) − 𝐴𝑤 ) 𝛿(𝐴𝑝 (𝑋) − 𝐴𝑝 )𝑒−𝑄𝑝 (𝑋)∕𝜃 , =𝛴. (87)
∫ 𝜕𝜃 𝐴 ,𝐴 𝑤 𝑝
𝑑𝑋
= 𝛿(𝑆𝑤 (𝑋) − 𝑆𝑤 ) 𝛿(𝜙(𝑋) − 𝜙)𝑒−𝑄𝑝 (𝑋)∕𝜃 , (80) We combine this equation with Eq. (59), taking the cross derivatives,
∫ 𝐴𝑝 𝐴
and finding
which we may write as ( ) ( )
𝜕𝛴 𝜕𝜇
= , (88)
1 𝜕𝐴𝑤 𝜃,𝐴 𝜕𝜃 𝐴 ,𝐴
𝑍(𝜃, 𝐴𝑤 , 𝐴𝑝 ) = 𝑒−𝑄𝐹 (𝜃,𝐴𝑤 ,𝐴𝑝 )∕𝜃 . (81) 𝑝 𝑤 𝑝
𝐴𝑝 𝐴
which may be written as
( ) ( )
We may repeat this procedure one more time. We rewrite the 𝜕𝜎 𝜕𝜇
partition function 𝑍(𝜃, 𝐴𝑤 , 𝐴𝑝 ), Eq. (80) as = . (89)
𝜕𝑆𝑤 𝜃 𝜕𝜃 𝑆
𝑤
+∞
𝑍(𝜃, 𝐴𝑤 , 𝐴𝑝 ) = 𝑑𝑄𝑝 𝑒−𝑄𝑝 ∕𝜃 𝑍(𝑄𝑝 , 𝐴𝑤 , 𝐴𝑝 ) , (82) 4. Fluctuations and agiture
∫−∞
where Our starting point are Eqs. (41) and (47), i.e.,

𝑍(𝑄𝑝 , 𝐴𝑤 , 𝐴𝑝 ) = 𝑑𝑋 𝑄𝐺 (𝜃, 𝜇, 𝜋) = −𝜃 ln[𝑍(𝜃, 𝜇, 𝜋)]

∫ [ ]
= −𝜃 ln 𝑑𝑋𝑒−𝑄𝑝 (𝑋)∕𝜃+𝜇𝐴𝑤 (𝑋)∕𝜃+𝜋𝐴𝑝 (𝑋)∕𝜃 . (90)
𝛿(𝑄𝑝 (𝑋) − 𝑄𝑝 ) 𝛿(𝐴𝑤 (𝑋) − 𝐴𝑤 ) 𝛿(𝐴𝑝 (𝑋) − 𝐴𝑝 ) . (83) ∫
This immediately gives
This microcanonical partition function (as 𝑄𝑝 is the control variable) is ( )
also the (unnormalized) density of states with respect to the variables 𝜕𝑄𝐺 (𝜃, 𝜇, 𝜋)
𝑄𝑝 , 𝐴𝑤 and 𝐴𝑝 . It demonstrates that all states 𝑋 with the same 𝑄𝑝 , 𝜕𝜇 𝜃,𝜋
𝐴𝑤 and 𝐴𝑝 are equally probable. This brings us back to our starting 1
= 𝑑𝑋𝐴𝑤 (𝑋)𝑒−𝑄𝑝 (𝑋)∕𝜃+𝜇𝐴𝑤 (𝑋)∕𝜃+𝜋𝐴𝑝 (𝑋)∕𝜃
point: the entropy has its maximum when all states that comply with 𝑍∫
the constraints (31), (32) and (33) are equally probable. = 𝐴𝑤 . (91)

8
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

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

𝑄𝑝 = 𝑄𝐴 𝐵
𝑝 + 𝑄𝑝 , (98)
𝐴 𝐵
𝛴 = 𝛴 +𝛴 , (99)
𝐴𝑤 = 𝐴𝐴 𝐵
𝑤 + 𝐴𝑤 , (100)
𝐴𝑝 = 𝐴𝐴 𝐵
𝑝 + 𝐴𝑝 . (101)
We write 𝑄𝑝 in terms of the control variables 𝛴, 𝐴𝑤 and 𝐴𝑝 ,
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 𝜃 𝐴 , 𝜇𝐴 , 𝜋 𝐴 and 𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 ) = 𝑄𝐴 𝐴 𝐴 𝐴 𝐵 𝐵 𝐵 𝐵
𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 ) + 𝑄𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 ) . (102)
∇𝑃 𝐴 , whereas in region B we have 𝜃 𝐵 , 𝜇 𝐵 , 𝜋 𝐵 and ∇𝑃 𝐵 . Steady state flow is attained
when 𝜃 𝐴 = 𝜃 𝐵 and 𝜇𝐴 = 𝜇𝐵 and 𝜋 𝐴 = 𝜋 𝐵 according to Eqs. (110), (112) and (112). Since the total volumetric flow rate is conserved in the plug, we must
have that fluctuations in it must be zero, i.e,
𝛿𝑄𝑝 = 𝛿𝑄𝐴 𝐵
𝑝 + 𝛿𝑄𝑝 = 0 . (103)
Taking the derivative a second time with respect to 𝜇 gives
( 2 ) The fluctuations in 𝑄𝐴 𝐵 𝐴
𝑝 and 𝑄𝑝 come from fluctuations in 𝛴 and
𝜕 𝑄𝐺 (𝜃, 𝜇, 𝜋) 𝛴 𝐵 , 𝐴𝐴 and 𝐴 𝐵 , and 𝐴𝐴 and 𝐴𝐵 . We express 𝛿𝑄𝐴 and 𝛿𝑄𝐵 in terms of
− 𝜃 𝑤 𝑤 𝑝 𝑝 𝑝 𝑝
𝜕𝜇2 𝜃,𝜋 the fluctuations of these variables,
1
= 𝑑𝑋𝐴2𝑤 (𝑋)𝑒−𝑄𝑝 (𝑋)∕𝜃+𝜇𝐴𝑤 (𝑋)∕𝜃+𝜋𝐴𝑝 (𝑋)∕𝜃 𝛿𝑄𝐴 𝐴 𝐴 𝐴
𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
𝑍∫ ( )
[ ]2 𝜕𝑄𝐴 𝐴 𝐴 𝐴
1 𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
− 𝑑𝑋𝐴𝑤 (𝑋)𝑒−𝑄𝑝 (𝑋)∕𝜃+𝜇𝐴𝑤 (𝑋)∕𝜃+𝜋𝐴𝑝 (𝑋)∕𝜃 = 𝛿𝛴 𝐴
𝑍∫ 𝜕𝛴 𝐴
𝐴𝐴 𝐴
𝑤 ,𝐴𝑝
( )
= ⟨𝐴2𝑤 ⟩ − 𝐴2𝑤 = 𝛥𝐴2𝑤 , (92) 𝜕𝑄𝐴 𝐴 𝐴 𝐴
𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
+ 𝛿𝐴𝐴
𝑤
where we have defined the fluctuations 𝛥𝐴2𝑤 . We now combine Eqs. (1) 𝜕𝐴𝐴
𝑤 𝛴 𝐴 ,𝐴𝐴
𝑝
and (7) to find ( )
𝜕𝑄𝐴 𝐴 𝐴 𝐴
𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
𝐴𝑝 𝐴𝑤 + 𝛿𝐴𝐴
𝑝 . (104)
𝐴𝑤 = 𝐴 = 𝐴𝜙𝑆𝑤 . (93) 𝜕𝐴𝐴
𝑝
𝐴 𝐴𝑝 𝛴 𝐴 ,𝐴𝐴
𝑤

Likewise, we have
This allows us to transform Eq. (92) into
( ) 𝛿𝑄𝐵 𝐵 𝐵 𝐵
𝜃 𝜕(𝜙𝑆𝑤 ) 𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
𝛥(𝜙𝑆𝑤 )2 = . (94) ( )
𝐴 𝜕𝜇 𝜃,𝜋 𝜕𝑄𝐵 𝐵 𝐵 𝐵
𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
= 𝛿𝛴 𝐵
We see that the agiture 𝜃 seems proportional to the fluctuations 𝜕𝛴 𝐵
𝐴𝐵 𝐵
𝑤 ,𝐴𝑝
𝛥(𝜙𝑆𝑤 )2 . However, due to the term (𝜕𝜙𝑆𝑤 ∕𝜕𝜇)𝜃,𝜋 , the relation between ( )
𝜕𝑄𝐵 𝐵 𝐵 𝐵
𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
them is complex. We may calculate the porosity fluctuations by taking + 𝛿𝐴𝐵
𝑤
the partition function 𝑍(𝜃, 𝜇, 𝜋), Eq. (75), as starting point. We find 𝜕𝐴𝐵
𝑤 𝛴 𝐵 ,𝐴𝐵
𝑝
( ) ( )
𝜃 𝜕𝜙 𝜕𝑄𝐵 𝐵 𝐵 𝐵
𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
𝛥𝜙2 = , (95) + 𝛿𝐴𝐵 (105)
𝐴 𝜕𝜋 𝜃,𝜇 𝑝 .
𝜕𝐴𝐵
𝑝 𝛴 𝐵 ,𝐴𝐵
𝑤
where
[ ] There is no production of cluster entropy as this entropy is at a
1 maximum. However, cluster entropy may move between regions A and
𝛥𝜙2 = ⟨𝐴𝑝 (𝑋)2 ⟩ − 𝐴2𝑝 . (96)
𝐴2 B. Hence, we have
We note that the porosity field 𝜙 is a property of the matrix and not
𝛿𝛴 = 𝛿𝛴 𝐴 + 𝛿𝛴 𝐵 = 0 . (106)
the flow. Hence, Eq. (95) gives a direct link between the agiture 𝜃 and
the flow pressure 𝜋, Furthermore, we keep the control variables 𝐴𝑤 and 𝐴𝑝 fixed. This
( ) is of course not possible experimentally. However, as a theoretical
𝜕𝜋
𝜃 = 𝐴𝛥𝜙2 . (97) concept, it is permissible. Hence, we have that
𝜕𝜙 𝜃,𝜇
𝛿𝐴𝑤 = 𝛿𝐴𝐴 𝐵
𝑤 + 𝛿𝐴𝑤 = 0 , (107)
5. Conditions for steady state in a heterogeneous plug 𝛿𝐴𝑝 = 𝛿𝐴𝐴
𝑝 + 𝛿𝐴𝐵
𝑝 =0. (108)
We now combine Eqs. (103)–(108) to find
We consider in the following the conditions for steady state in a
heterogeneous plug. 𝛿𝑄𝑝 (𝛴, 𝐴𝑤 , 𝐴𝑝 )
( ) ( 𝐵)
We show in Fig. 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 𝑄𝑝 , 𝐴𝑤 , 𝐴𝑝 ⎢ 𝜕𝐴𝑤 𝐴 𝜕𝐴𝑤𝐵 ⎥ 𝑤
⎣ 𝛴 𝐴 ,𝐴𝐴
𝑝 𝑝 ⎦
𝛴 𝐵 ,𝐴𝐵
and 𝛴 that describe the flow in the plug. Strictly speaking, we should ( ) ( 𝐵)
⎡ 𝜕𝑄𝐴 𝜕𝑄𝑝 ⎤
change our notation, e.g., 𝑄𝑝 → 𝑝 , since we are considering the + ⎢ 𝑝
− ⎥ 𝛿𝐴𝐴 = 0 . (109)
entire plug. However, in order to avoid complicating the notation and ⎢ 𝜕𝐴𝐴 𝜕𝐴𝐵 ⎥ 𝑝
⎣ 𝑝 𝛴 𝐴 ,𝐴𝐴
𝑤
𝑝 𝛴 𝐵 ,𝐴𝐵 ⎦
𝑤

9
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

We now assume that the cluster entropy fluctuations 𝛿𝛴 𝐴 are indepen- statistical mechanics, and local steady state flow plays the role of
dent of the fluctuations in 𝛿𝐴𝐴 𝐴
𝑤 and 𝛿𝐴𝑝 . This leads to 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. (2018), rendering it a com-
𝜕𝛴 𝐴 plete pseudo-thermodynamics in the spirit of Edwards and Oakeshot’s
𝐴𝐴 𝐴
𝑤 ,𝐴𝑝
( 𝐵 𝐵 𝐵 𝐵 ) pseudo-thermodynamics for powders (Edwards and Oakeshott, 1989).
𝜕𝑄𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 )
= = 𝜃𝐵 , (110) 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 hydro-
where we have used Eqs. (49) and (55). dynamic description to a continuum scale physics, thus solving the
If we now furthermore assume that the fluctuations in 𝛿𝐴𝐴
𝑤 and 𝛿𝐴𝑝
𝐴
up-scaling problem. The key link is the partition function 𝑍(𝜃, 𝜇, 𝜋)
are uncorrelated, we find defined in Eq. (41). The integral is over all physical fluid configurations,
( 𝐴 𝐴 𝐴 𝐴 )
𝜕𝑄𝑝 (𝛴 , 𝐴𝑤 , 𝐴𝑝 ) and the physics sits in the measure 𝑑𝑋 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
= = 𝜇𝐵 , (111) variables admits an immense level of ignorance. This ignorance is
𝜕𝐴𝐵𝑤 𝐵 𝐵 𝛴 ,𝐴𝑝
reflected in the pseudo-thermodynamics as the cluster entropy. We call
and 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
= 𝐵
= 𝜋𝐵 . (112) 𝜃 may be measured, or controlled, experimentally or computationally.
𝜕𝐴𝑝
𝛴 𝐵 ,𝐴𝐵
𝑤 The flow pressure 𝜋 seems also difficult to measure. Measuring the flow
These three criteria for steady state flow, Eqs. (110), (111) and derivative is, on the other hand, seemingly easier to measure. In terms
(112), are analogous to the criteria for two open systems in thermal of the average seepage velocity and the saturation, it may be written
contact and in equilibrium: here the temperature, pressure and the ( )
𝜕𝑣𝑝 (𝜎, 𝑆𝑤 )
chemical potential must be the same. 𝜇= . (113)
𝜕𝑆𝑤 𝜎
It is important to point out the following. In ordinary thermodynam-
ics, the rule is that the conjugate of a conserved quantity is constant in a Since we are assuming steady-state flow, the cluster entropy is at a
heterogeneous system at equilibrium. This is e.g., the argument for the maximum and therefore constant.
temperature being the same everywhere in the system at equilibrium We also note that we have used five different ensembles in this
as the entropy is at a maximum. This is the same argument as we have work, where control parameters are respectively
used here which leads to Eq. (110). However, consider two magnets (𝜃, 𝜇, 𝜋) ,
with different magnetic susceptibilities in contact which is placed in
(𝜃, 𝜇, 𝐴𝑝 ) ,
a uniform magnetic field 𝐻. The conjugate of the magnetization 𝑀,
which is the magnetic flux 𝐵, is not equal in the two magnets in (𝛴, 𝜇, 𝐴𝑝 ) ,
contact, 𝐵 𝐴 ≠ 𝐵 𝐵 . 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., 𝛿𝑀 𝐴 = −𝛿𝑀 𝐵 while simultaneously keeping the total
internal energy and the total entropy constant, i.e., 𝛿𝑈 = 𝛿𝑈 𝐴 +𝛿𝑈 𝐵 = 0 The first of these ensembles is easily realizable experimentally as this
and 𝛿𝑆 = 𝛿𝑆 𝐴 + 𝛿𝑆 𝐵 = 0. In our system — steady-state two-phase flow describes an REA communicating freely with the rest of the porous
in a porous medium — the following question then becomes central: is medium. The second one, with 𝐴𝑝 controlled, is feasible e.g. by using
it possible to fulfill 𝛿𝑄𝑝 = 0 and 𝛿𝛴 = 0, and at the same time have 3D printing techniques to construct the porous medium. The two
𝛿𝐴𝐴 𝐵 𝐴 𝐵
𝑤 = −𝛿𝐴𝑤 and 𝛿𝐴𝑝 = −𝛿𝐴𝑝 ? The answer to this question is not
last three ensembles must be seen as theoretical constructs. We note,
obvious. Only experiments or computations on models may provide however, that this is standard procedure in thermodynamics in that one
an answer. If the answer is no, only Eq. (110) will be valid, and not always starts with the Gibbs relation, which relates variations in the
Eqs. (111) and (112). internal energy to variations in the extensive variables irrespective of
whether this ensemble is accessible experimentally or not.
6. Discussion and conclusion We see that a steady-state condition in a heterogeneous porous
medium such as that shown in Fig. 4 is that the agiture is con-
A central problem in the physics of porous media is how to scale stant throughout the entire system, Eq. (110). On the other hand, the
up knowledge of the physics at the pore level to the continuum scale, two additional equilibrium conditions, Eqs. (111) and (112), hinge
often referred to as the Darcy scale, where the pores are negligible on additional assumptions that need to be verified experimentally or
in size. The notion of up-scaling is of course only meaningful if we computationally. In light of the above discussion, verifying the validity
know which coarse-scale physics we are up-scaling to, and the tradi- of Eq. (111) is the easier one of the three equations.
tional relative permeability theory has well known weak points. The The probability to find a flow configuration 𝑝(𝑋) in the REA which
theory of Hansen et al. (2018) is a radically different approach to the is the central to the theory we present here. It was emphasized at the
coarse-scale physics. This theory is formally similar to some parts of beginning of Section 3.3 that the fluid configuration 𝑋 refers to the
thermodynamics, but it is the flow rates, and the relations between flow configuration in the REA, not the entire plane cutting through the plug,
rates, that are the central players. In its core the theory is based on the  = 𝑋 ∪ 𝑋𝑟 . It is crucial for the theory that 𝑝(𝑋) does not depend on
volumetric flow rate being an Euler homogeneous function. the properties of the entire plug. Recent numerical work by Fyhn et al.
In the present paper, we have developed a statistical mechan- indicates that this is indeed correct.
ics framework for immiscible incompressible two-phase flow. In this Given all these caveats and open question, the fact remains that the
framework, total flow rate plays the same role as energy in ordinary framework we have presented here, opens up for viewing immiscible

10
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

two-phase flow in porous media in a different way than before. We have Gao, Y., Lin, Q., Bijeljic, B., Blunt, M.J., 2020. Pore-scale dynamics and the multiphase
here divided the problem into 1. constructing a framework by which Darcy law. Phys. Rev. Fluids 5, 013801. http://dx.doi.org/10.1103/PhysRevFluids.
5.013801.
the concepts that are necessary are put in place, and 2. connecting these Gray, W.G., Miller, C.T., 2014. Introduction to the Thermodynamically Constrained
concepts to the underlying fluid dynamics and thermodynamics. This Averaging Theory for Porous Medium Systems. Springer Verlag, Berlin, http://dx.
paper accomplishes the first goal. doi.org/10.1007/978-3-319-04010-3.
Grøva, M., Hansen, A., 2011. Two-phase flow in porous media: power-law scaling of
effective permeability. J. Phys.: Conf. Ser 319, 012009. http://dx.doi.org/10.1088/
Declaration of competing interest 1742-6596/319/1/012009.
Hansen, A., Ravndal, F., 1981. Klein’s paradox and its resolution. Phys. Scr. 23, 1036.
The authors declare the following financial interests/personal re- http://dx.doi.org/10.1088/0031-8949/23/6/002.
Hansen, A., Sinha, S., Bedeaux, D., Kjelstrup, S., Gjennestad, M.A., Vassvik, M., 2018.
lationships which may be considered as potential competing inter-
Relations between seepage velocities in immiscible, incompressible two-phase flow
ests: Alex Hansen reports financial support was provided by Research in porous media. Transp. Por. Med 125, 565. http://dx.doi.org/10.1007/s11242-
Council of Norway. Alex Hansen reports a relationship with Research 018-1139-6.
Council of Norway that includes: funding grants. Hassanizadeh, S.M., Gray, W.G., 1990. Mechanics and thermodynamics of multiphase
flow in porous media including interphase boundaries. Adv. Water Res 13, 169.
http://dx.doi.org/10.1016/0309-1708(90)90040-B.
Data availability Hassanizadeh, S.M., Gray, W.G., 1993a. Towards an improved description of the
physics of two-phase flow. Adv. Water Res 16, 53. http://dx.doi.org/10.1016/0309-
No data was used for the research described in the article. 1708(93)90029-F.
Hassanizadeh, S.M., Gray, W.G., 1993b. Thermodynamic basis of capillary pres-
sure in porous media. Water Resour. Res 29, 3389. http://dx.doi.org/10.1029/
Acknowledgments 93WR01495.
Hilfer, R., 2006. Capillary pressure, hysteresis and residual saturation in porous media.
The authors thank Dick Bedeaux, Carl Fredrik Berg, Daan Frenkel, Physica A 359, 119. http://dx.doi.org/10.1016/j.physa.2005.05.086.
Hilfer, R., 2006a. Macroscopic capillarity and hysteresis for flow in porous media. Phys.
Hursanay Fyhn, Signe Kjelstrup, Marcel Moura, Knut Jørgen Måløy, Rev. E 73, 016307. http://dx.doi.org/10.1103/PhysRevE.73.016307.
Håkon Pedersen, Subhadeep Roy, Ole Torsæter and Øivind Wilhelmsen Hilfer, R., 2006b. Macroscopic capillarity without a constitutive capillary pressure
for interesting discussions. This work was partly supported by the function. Physica A 371, 209. http://dx.doi.org/10.1016/j.physa.2006.04.051.
Research Council of Norway through its Center of Excellence funding Hilfer, R., Besserer, H., 2000. Macroscopic two-phase flow in porous media. Physica B
279, 125. http://dx.doi.org/10.1016/S0921-4526(99)00694-8.
scheme, project number 262644. Hilfer, R., Döster, F., 2010. Percolation as a basic concept for capillarity. Transp. Por.
Med 82, 507. http://dx.doi.org/10.1007/s11242-009-9395-0.
References Jaynes, E.T., 1957. Information theory of statistical mechanics. Phys. Rev 106, 620.
http://dx.doi.org/10.1103/PhysRev.106.620.
Kjelstrup, S., Bedeaux, D., Hansen, A., Hafskjold, B., Galteland, O., 2018. Non-
Armstrong, R.T., McClure, J.E., Robins, V., Liu, Z., Arns, C.H., Schlüter, S., Berg, S.,
isothermal transport of multi-phase fluids in porous media. The entropy production.
2019. Porous media characterization using Minkowski functionals: theories, appli-
Front. Phys 6, 126. http://dx.doi.org/10.3389/fphy.2018.00126.
cations and future directions. Transp. Porous Media 130, 305. http://dx.doi.org/
Kjelstrup, S., Bedeaux, D., Hansen, A., Hafskjold, B., Galteland, O., 2019. Non-
10.1007/s11242-018-1201-4.
isothermal transport of multi-phase fluids in porous media. Constitutive equations.
Aursjø, O., Erpelding, M., Tallakstad, K.T., Flekkøy, E.G., Hansen, A., Måløy, K.J.,
Front. Phys 6, 150. http://dx.doi.org/10.3389/fphy.2018.00150.
2014. Film flow dominated simultaneous flow of two viscous incompressible fluids
Knudsen, H.A., Hansen, A., 2006. Two-phase flow in porous media: dynamical phase
through a porous medium. Front. Phys 2, 63. http://dx.doi.org/10.3389/fphy.2014.
transition. Eur. Phys. J. B 49, 109. http://dx.doi.org/10.1140/epjb/e2006-00019-y.
00063.
Lanza, F., Hansen, A., Rosso, A., Talon, L., 2022. Transp. Por. Med 145, 245. http:
Barenblatt, G.I., Patzek, T.W., Silin, B.B., 2002. The mathematical model of non-
//dx.doi.org/10.1007/s11242-022-01848-7,
equilibrium effects in water-oil displacement. SPE-75169-MS http://dx.doi.org/10.
Måløy, K.J., Moura, M., Hansen, A., Flekkøy, E.G., Toussaint, R., 2021. Burst dynamics,
2118/75169-MS.
upscaling and dissipation of slow drainage in porous media. Front. Phys 9, 796019.
Bear, J., 1988. Dynamics of Fluids in Porous Media. Dover, Mineola, http://dx.doi.org/
http://dx.doi.org/10.3389/fphy.2021.796019.
10.1097/00010694-197508000-00022. McClure, J.E., Armstrong, R.T., Berg, S., Geometric evolution as a source of discontin-
Bear, J., Bachmat, Y., 2012. Introduction to Modeling of Transport Phenomena in uous behavior in soft condensed matter, arXiv:1906.04073; http://dx.doi.org/10.
Porous Media. Springer, Berlin, http://dx.doi.org/10.1007/978-94-009-1926-6. 48550/arXiv.1906.04073.
Bedeaux, D., Kjelstrup, S., 2022. Fluctuation-dissipiation theorems for multiphase flow McClure, J.E., Armstrong, R.T., Berrill, M.A., Schlüter, S., Berg, S., Gray, W.G.,
in porous media. Entropy 24, 46. http://dx.doi.org/10.3390/e24010046. Miller, C.T., 2018. Geometric state function for two-fluid flow in porous media.
Blunt, M.J., 2017. Multiphase Flow in Permeable Media. Cambridge Univ. Press, Phys. Rev. Fluids 3, 084306. http://dx.doi.org/10.1103/PhysRevFluids.3.08430.
Cambridge, http://dx.doi.org/10.1017/9781316145098. McClure, J.E., Berg, S., Armstrong, R.T., 2021a. Capillary fluctuations and energy
Callen, H.B., 1974. Thermodynamics as a science of symmetry. Found. Phys 4, 423. dynamics for flow in porous media. Phys. Fluids 33, 083323. http://dx.doi.org/
http://dx.doi.org/10.1007/BF00708519. 10.1063/5.0057428.
Callen, H.B., 1991. Thermodynamics and an Introduction To Thermostatistics, 2nd McClure, J.E., Berg, S., Armstrong, R.T., 2021b. Thermodynamics of fluctuations based
Edition Wiley, New York, ISBN: 978-0-471-86256-7. on time-and-space averages. Phys. Rev. E 104, 035106. http://dx.doi.org/10.1103/
Cheon, H.L., Fyhn, H., Hansen, A., Wilhelmsen, Ø., Sinha, S., Steady-state two-phase PhysRevE.104.035106.
flow of compressible and incompressible fluids in a capillary tube of varying radius, Niessner, J., Berg, S., Hassanizadeh, S.M., 2011. Comparison of two-phase Darcy’s
arXiv:2207.10503; https://doi.org/10.48550/arXiv.2207.10503. law with a thermodynamically consistent approach. Transp. Por. Med 88, 133.
Döster, F., Hönig, O., Hilfer, R., 2012. Horizontal flow and capillarity-driven redis- http://dx.doi.org/10.1007/s11242-011-9730-0.
tribution in porous media. Phys. Rev. E 86, 016317. http://dx.doi.org/10.1103/ Roy, S., Pedersen, H., Sinha, S., Hansen, A., 2022. The co-moving velocity in immiscible
PhysRevE.86.016317. two-phase flow in porous media. Transp. Por. Med 143, 69. http://dx.doi.org/10.
Edwards, S.F., Oakeshott, R.B.S., 1989. Theory of powders. Physica A 157, 1080. 1007/s11242-022-01783-7.
http://dx.doi.org/10.1016/0378-4371(89)90034-4. Roy, S., Sinha, S., Hansen, A., Immiscible two-phase flow in porous media: effective
Erpelding, M., Sinha, S., Tallakstad, K.T., Hansen, A., Flekkøy, E.G., Måløy, K.J., 2013. rheology in the continuum limit, arXiv:1912.05248; https://doi.org/10.48550/
History independence of steady state in simultaneous two-phase flow through two- arXiv.1912.05248.
dimensional porous media. Phys. Rev. E 88, 053004. http://dx.doi.org/10.1103/ Roy, S., Sinha, S., Hansen, A., 2019. Effective rheology of two-phase flow in a capillary
PhysRevE.88.053004. fiber bundle model. Front. Phys 7, 92. http://dx.doi.org/10.3389/fphy.2019.00092.
Feder, J., Flekkøy, E.G., Hansen, A., 2022. Physics of Flow in Porous Media. Cambridge Roy, S., Sinha, S., Hansen, A., 2020. Flow-area relations in immiscible two-phase flow
Univ. Press, Cambridge, http://dx.doi.org/10.1017/9781009100717. in porous media. Front. Phys 8, 4. http://dx.doi.org/10.3389/fphy.2020.00004.
Feynman, R.P., 1948. A relativistic cut-off for classical electrodynamics. Phys. Rev 74, Sahimi, M., 2011. Flow and Transport in Porous Media and Fractured Rock: From
939. http://dx.doi.org/10.1103/PhysRev.74.939. Classical Methods to Modern Approaches. Wiley, New York, http://dx.doi.org/10.
Fyhn, H., Sinha, S., Hansen, A., Local statistics of immiscible and incompressible two- 1002/9783527636693.
phase flow in porous media, arXiv:2209.00030; http://dx.doi.org/10.48550/arXiv. Shannon, C.E., 1948. A mathematical theory of communication. Bell Syst. Tech. J. 27,
2209.00030. 379. http://dx.doi.org/10.1002/j.1538-7305.1948.tb01338.x.
Fyhn, H., Sinha, S., Roy, S., Hansen, A., 2021. Rheology of immiscible two-phase flow in Sinha, S., Bender, A.T., Danczyk, M., Keepseagle, K., Prather, C.A., Bray, J.M.,
mixed wet porous media: dynamic pore network model and capillary fiber bundle Thrane, L.W., Seymour, J.D., Codd, S.L., Hansen, A., 2017. Effective rheology of
model results. Transp. Por. Med 139, 491. http://dx.doi.org/10.1007/s11242-021- two-phase flow in three-dimensional porous media: experiment and simulation.
01674-3. Transp. Por. Med 119, 77–94. http://dx.doi.org/10.1007/s11242-017-0874-4.

11
A. Hansen et al. Advances in Water Resources 171 (2023) 104336

Sinha, S., Hansen, A., 2012. Effective rheology of immiscible two-phase flow in Valavanides, M.S., Constantinides, G.N., Payatakes, A.C., 1998. Mechanistic model of
porous media. Europhys. Lett. 99, 44004. http://dx.doi.org/10.1209/0295-5075/ steady-state two-phase flow in porous media based on ganglion dynamics. Transp.
99/44004. Por. Med 30, 267. http://dx.doi.org/10.1023/A:1006558121674.
Sinha, S., Hansen, A., Bedeaux, D., Kjelstrup, S., 2013. Effective rheology of bubbles Wang, Y., Aryana, S.A., Allen, M.B., 2019. An extension of Darcy’s law incorporating
moving in a capillary tube. Phys. Rev. E 87, 025001. http://dx.doi.org/10.1103/ dynamic length scales. Adv. Water Res. 129, 70. http://dx.doi.org/10.1016/j.
PhysRevE.87.025001. advwatres.2019.05.010.
Tallakstad, K.T., Knudsen, H.A., Ramstad, T., Løvoll, G., Måløy, K.J., Toussaint, R., Whitaker, S., 1986. Flow in porous media II: The governing equations for immiscible,
Flekkøy, E.G., 2009a. Steady-state two-phase flow in porous media: statistics and two-phase flow. Transp. Por. Med 1, 105. http://dx.doi.org/10.1007/BF00714688.
transport properties. Phys. Rev. Lett. 102, 074502. http://dx.doi.org/10.1103/ Wyckoff, R.D., Botset, H.G., 1936. The flow of gas-liquid mixtures through
PhysRevLett.102.074502. unconsolidated sands. Physics 7, 325. http://dx.doi.org/10.1063/1.1745402.
Tallakstad, K.T., Løvoll, G., Knudsen, H.A., Ramstad, T., Flekkøy, E.G., Måløy, K.J., Xu, X., Wang, X., 2014. Non-Darcy behavior of two-phase channel flow. Phys. Rev. E
2009b. Steady-state simultaneous two-phase flow in porous media: an experimental 90, 023010. http://dx.doi.org/10.1103/PhysRevE.90.023010.
study. Phys. Rev. E 80, 036308. http://dx.doi.org/10.1103/PhysRevE.80.036308. Yiotis, A.G., Dollari, A., Kainourgiakis, M.E., Salin, D., Talon, L., 2019. Nonlinear Darcy
Valavanides, M.S., 2012. Steady-state two-phase flow in porous media: review of flow dynamics during ganglia stranding and mobilization in heterogeneous porous
progress in the development of the DeProF theory bridging pore- to statistical domains. Phys. Rev. Fluids 4, 114302. http://dx.doi.org/10.1103/PhysRevFluids.4.
thermodynamics-scales. Oil Gas Sci. Technol 67, 787. http://dx.doi.org/10.2516/ 114302.
ogst/2012056. Zhang, Y., Bijeljic, B., Gao, Y., Lin, Q., Blunt, M.J., 2021. Quantification of nonlinear
Valavanides, M.S., 2018. Review of steady-state two-phase flow in porous media: multiphase flow in porous media. Geophys. Res. Lett 48, e2020GL090477. http:
independent variables, universal energy efficiency map, critical flow conditions, //dx.doi.org/10.1029/2020GL090477.
effective characterization of flow and pore network. Transp. Porous Media 123,
45. http://dx.doi.org/10.1007/s11242-018-1026-1.

12