Capillary energy barriers have important consequences for immiscible fluid flow in porous media. We derive a time-and-space averaging theory to account for the non-equilibrium behavior and understand the role of athermal capillary fluctuations in the context of their relationship to larger scale phenomenological equations. The formulation resolves several key challenges associated with two-fluid flow in porous media: (1) geometric and thermodynamic quantities are constructed as smooth functions of time based on time-and-space averages; (2) averaged thermodynamics are developed for films; (3) multi-scale fluctuation terms are identified, which account for transient behaviors of interfaces and films that occur due to pore-scale events; (4) geometric constraints are derived and imposed on the averaged thermodynamics; (5) a new constitutive model is proposed for capillary pressure dynamics that includes contributions from films; and (6) a time-and-space criterion for representative elementary volume is established based on capillary fluctuations. Capillary fluctuations are assessed quantitatively based on pore-scale simulations and experimental core-flooding data.

## I. INTRODUCTION

The non-equilibrium response to fluid flow through porous materials is inextricably linked to the length and time scales for processes that occur within the solid micro-structure.^{1–3} It is well known that fluid pressures fluctuate during immiscible displacement as a consequence of pore-scale events.^{4–6,72} The associated fluctuations are athermal and cooperative in nature based on the fact that they arise due to the influence of capillary forces.^{7} Theoretical treatment of these fluctuations is of central importance to modeling due to the fact that athermal fluctuations may not obey detailed balance.^{8} Since the reversibility of molecular interactions is central to the development of Onsager's non-equilibrium theory, the validity of phenomenological equations depends on the properties of fluctuations within the system.^{9,10} The validity of near-equilibrium approximations and Onsager's theory has been an important assumption in many theoretical developments for fluid flow in porous media.^{11–15}

When multiple length scales are present within a system, accompanying timescales will arise when considering the system dynamics. The linkage between temporal, spatial, and energy scales is summarized in Fig. 1. Since information cannot propagate instantaneously, processes that occur over small length scales will tend to be more rapid compared to those that occur over larger length scales. Given a particular length scale, the energy scale relates to the work needed to move matter and energy at the requisite distance. The rate that information propagates in the system is embedded in the phenomenological coefficients associated with particular physical processes. With respect to the validity of linear response theory, Onsager's theory holds remarkably well in small systems, implying that local near-equilibrium approximations are usually reasonable.^{16} At small length scales, near-equilibrium approximations are expected to fail only in extreme cases, such as in the immediate vicinity of shock waves.^{17} However, as the length scale associated with a process increases, so too does the required timescale to reach equilibrium. For immiscible fluids in a geological setting, the timescale required for diffusive processes to reach equilibrium is remarkably slow; the timescale to reach equilibrium for Ostwald ripening in CO_{2}–brine systems can be years.^{18} Within this context, slow physics represent a critical consideration when applying near-equilibrium assumptions at macroscopic length scales. The slow relaxation of diffusive systems explains the near-ubiquitous failure of Fick's law in porous media.^{19–22} While the timescale associated with capillary fluctuations is measured in seconds, the timescale to reach equilibrium can be hours to days in experimental systems.^{23} It is therefore of great interest to develop non-equilibrium theory that does not rely on macroscopic near-equilibrium assumptions. In this paper, we develop such a theory.

For multiphase flows in porous media, we identify three length scales of primary interest: (1) the interfacial length scale associated with the width of the interface (10–50 nm); (2) the pore length scale associated with the solid micro-structure (∼mm); and (3) the Darcy length scale associated with a macroscopic flow (cm to km). For each length scale, particular processes of interest can be identified. At the interface length scale, film swelling effects can occur due to flow; coalescence events occur based on the merging of fluid regions, which form loops or otherwise alter the fluid topology.^{24} Fluid singularities result from the associated disruption to the interface mean curvature.^{25–38} Snap-off events are also common in typical flow processes.^{39} At the pore length scale, Haines jumps occur when fluid spontaneously fills regions of the pore structure as fluids migrate between capillary energy barriers.^{40} These microscopic events are a primary source of non-equilibrium behavior and are therefore a focus area when considering the physical meaning of near equilibrium approximations. The length scale of interest for flow processes is the Darcy scale, typically orders of magnitude larger than that of the pore scale.

In Sec. I A, we provide a conceptual overview of the basic physical considerations needed to model immiscible displacement in porous media. Formal theoretical development proceeds in six steps: (1) definition of the thermodynamic system; (2) averaging in time-and-space; (3) rate of external work in Darcy-scale systems; (4) analysis for fluctuations in solid wetting energy; (5) derivation of geometric constraints; and (6) formulation of flux–force entropy and derivation of phenomenological equations for two-phase flows. The derivation leads to explicit criteria for the validity of phenomenological equations based on the structure of capillary fluctuations in the system. Using results from simulation and experimental data, we directly evaluate the derived capillary fluctuation terms and consider the consequences for definition of representative elementary volume (REV) and the development of averaged descriptions for two-fluid flows in porous media.

### A. Background

At the pore scale, capillary forces are dominant and determine how fluids are distributed within the solid micro-structure. Capillary pressure accounts for the contribution of surface forces created by the surface energy and meniscus curvature.^{41} The Young–Laplace equation relates the pressure difference between the adjoining fluids and the capillary pressure given by the product of the interfacial tension *γ* and mean curvature,

where *R*_{1} and *R*_{2} are the principal curvature radii at points on the fluid meniscus. The Laplace equation holds only at mechanical equilibrium and does not account for dynamic processes that involve moving interfaces. This necessarily includes any process where the fluid volumes change, since these changes can only occur based on movement of the interface. Situations also arise where interface movement occurs at constant volume, such as capillary waves. To describe these cases, it is necessary to understand the non-equilibrium response of the fluid pressures based on the capillary dynamics within the system.

As fluids squeeze through the micro-structure of complex materials subject to the influence of capillary forces, rapid changes to the fluid configuration occur as the system jumps between metastable states. The displacement can be sub-divided into a sequence of *isons* and *rheons.*^{4,5,43} The associated events are depicted in Fig. 2. *Isons* correspond to reversible displacement events where the fluid meniscus invades a pore throat under quasi-static conditions, i.e., the fluid pressure difference is balanced by the capillary forces in the throat. The fluid pressure difference increases proportionately to the meniscus curvature as the interface is pushed into increasingly narrow parts of the pore throat. When the meniscus moves through the narrowest part of the pore neck, pore-filling occurs spontaneously based on an event often known as a Haines jump.^{40} As the meniscus curvature decreases, the fluid pressure difference significantly exceeds the capillary pressure based on the interface curvature and fluid rapidly fills the pore. This irreversible event is known as a *rheon*. Pore-filling events that transpire as Haines jumps can dissipate a significant amount of energy.^{1,2} Furthermore, they occur very quickly (milliseconds to seconds) compared to typical laboratory- and reservoir-scale process (hours to weeks). Considering large systems, processes of interest involve very large numbers of such events. Understanding their overall effect is of central importance to understanding transient effects for multiphase flows in porous media.

The basis to consider multiphase flow in porous media as a non-ergodic phenomenon is due to the rapid rate for local events at the pore-scale. In ergodic systems, energetic micro-states are considered to be locally well-mixed due to thermal effects. Since the operative mixing mechanism is diffusive, the length scale for thermal mixing can be estimated based on the self-diffusion coefficient and timescale.^{50} When a mechanical event transfers energy more rapidly than diffusive thermal mechanisms, local symmetry breaking may occur. Figure 3 considers the length and time scale for seven Haines jump events observed during displacement in a microfluidic system.^{1} The distance traveled by the fluid meniscus is plotted vs the elapsed time for each event. This can be compared to the length scale for diffusive mixing during the elapsed time is $\Delta t$, given by

where $D=2.299\xd710\u22125$ cm^{2}/s is the self-diffusion coefficient for water. We classify a phenomenon as ergodic if the distance traveled over a particular time interval is smaller than the diffusive length scale. A phenomenon is classified as non-ergodic if the distance traveled during an interval of time exceeds this distance. Figure 3 clearly shows that Haines jumps begin in the ergodic region but cross into the non-ergodic regime as the capillary energy is converted into kinetic energy. Over the timescale for a Haines jump, the length scale for thermal mixing is an order of magnitude smaller than the characteristic length scale for the event. The system can therefore be considered to be thermally well-mixed at the scale of $D\Delta t$, but not at the scale of a Haines jump.

At the Darcy scale, energy that is dissipated from pore-scale events can be estimated from basic concepts of work and energy. A typical experiment can be treated as a system where external work is performed in two primary ways. First is the work associated with fluid flow. The volumetric flow rate is $Qi$ for each fluid $i\u2208{w,n}$ results due to external work performed on the system by some external potential $\Phi i$, e.g., due to an external pressure gradient or body force. Since work is defined as the applied force multiplied by the distance, the total rate of work is

where $\Delta \Phi i=\Phi i(in)\u2212\Phi i(out)$ is the potential difference between the inlet and the outlet and $nb$ is the boundary normal vector. This expression is directly implied from the definition of mechanical work. The boundary terms can also be obtained by applying the divergence theorem in consideration of the internal stresses.^{43,65} Conceptually, the work due to Darcy flow occurs in a steady-state system, meaning that there are no net changes to the internal energy, with each fluid forming static connected pathways through the sample that allow flow to occur. The conventional relative permeability relationship predicts the flux per area as

where $qi=Qi/A$, *k _{ri}* is the relative permeability,

*μ*is the dynamic viscosity, and

_{i}**K**is the permeability tensor for the material. This can be combined with Eq. (3) to approximate the rate of work based on an applied external potential gradient. The rate of energy dissipated within the system is thereby embedded in

*k*.

_{ri}Pressure–volume work occurs when one fluid displaces another. Considering a wetting and non-wetting fluid with pressure *p _{w}* and

*p*, this contribution to the external work is given by

_{n}External work performed on the system is transferred to the internal energy modes, which include contributions to the internal thermal energy from dissipated heat. Changes to the surface energy account for a significant part of the external pressure–volume work. If the solid material is water–wet, the change in surface energy is proportional to the change in surface area based on the interfacial tension *γ _{wn}*. Morrow, therefore, defined the displacement efficiency as the fraction of pressure–volume work that is converted to surface energy.

^{5}Seth and Morrow showed that the efficiency for a primary drainage varies from 10% to 95% depending on the material type.

^{42}The dissipated energy can also be estimated from the mechanistic picture illustrated in Fig. 2. The pressure–volume work associated with an initial

*ison*can be considered to be reversible based on the fact that the displacement is quasi-static. Since a

*rheon*is spontaneous, the work associated with the ensuing sub-ison accounts for the dissipated energy after any surface energy changes generated by the event has been subtracted. With energy input into the system given by Eqs. (4) and (5), the internal energy dynamics of the system can be treated based on the conservation of energy.

Models for capillary pressure dynamics have been previously developed to account for the macroscopic behavior during displacement. The first such models were developed using concepts of work and energy.^{5} Subsequent models incorporated concepts of non-equilibrium thermodynamics^{11,15} and kinematics.^{13,44} Previously developed models fail to deal with the challenges associated with the non-smooth response of the fluid pressure and interfacial curvature. This critical challenge is associated with changes to the fluid topology, which necessarily cause a discontinuity in the interfacial curvature.^{24,30,37} The apparent geometric discontinuity links to questions about whether or not multi-fluid systems can be considered to be in a near-equilibrium state. Bear and Nitao argue that while relatively simple macroscopic thermodynamic equations of state must exist at equilibrium due to the Gibbs phase rule, these relations will break down under non-equilibrium conditions if the mean square deviation is too large.^{14} For systems with slow physics, near equilibrium conditions can break down as the length scale for the REV increases due to the larger timescale required for mixing to occur. We argue that the breakdown of near-equilibrium behavior will occur at large length scales based on the argument summarized in Fig. 3. Near-equilibrium conditions should hold at spatial scales smaller than $D\Delta t$, with $\Delta t$ being the timescale for pore-scale events. Averaging in time-and-space provides a way to upscale the formulation to larger length and time scales based on a near-equilibrium approximation at a small scale. Upscaling is then performed based on integrals, which do not require any near-equilibrium approximation at the larger scale. In this work, we consider how multi-scale effects influence the transient capillary pressure behavior and consider consequences for the rate of energy dissipation. Multi-scale fluctuation terms are identified explicitly, which can be directly evaluated and compared to other energy dynamics within the system. Conditions derived based on these fluctuation terms clearly demonstrate when Onsager theory can be applied to derive phenomenological models such as the conventional relative permeability relationship. A new expression for the dynamic capillary pressure is also derived, which directly includes effects due to fluctuations. Capillary fluctuations are then directly evaluated based on pore-scale simulations and experimental data.

## II. THEORETICAL APPROACH

In this section, we present the theoretical approach used to derive constitutive models for immiscible fluid flow through porous media. Our approach relies on averaging in both time and space, which is distinct from previously developed models. Capillary fluctuation terms are included explicitly, and special treatment is needed to develop a flux–force entropy inequality that accounts for the associated energy dynamics. The flux–force form of the entropy inequality can be exploited to derive standard expressions for relative permeability as well as a dynamic capillary pressure relationship to describe unsteady displacement. The theoretical approach is organized into six steps:

**Thermodynamic description**—variables required to formulate a classical thermodynamic description at a scale where the behavior is ergodic, including effects due to fluid phases, interfaces, and films;**Averages in space and time**—averaging is applied to develop thermodynamic forms that hold at larger length and time scales. Fluctuation terms arise as a consequence of averaging in both space and time.**Rate of work and entropy inequality**—an entropy inequality is developed by defining the rate of external work associated with fluid flow;**Surface wetting fluctuations**—analysis is performed to interpret the surface wetting fluctuation term;**Geometric constraints**—a geometric constraint is needed to relate the time derivative for the interfacial area to the time derivative for the fluid volume fractional. A relationship between partial derivatives for geometric quantities is derived;**Flux–force form of entropy inequality**—all results are combined to form a flux–force form of the entropy inequality. A REV constraint is derived for the fluctuations. Provided the constraint is met, standard relative permeability relationships are obtained as well as a form to approximate the capillary pressure dynamics during unsteady displacement.

By explicitly considering the contribution from fluctuations in the flux–force form of the entropy inequality, conditions for the validity of Darcy-scale phenomenological equations can be developed explicitly. The derived condition is a constraint on the REV for the system, which considers both spatial and temporal aspects. The theory does not rely on the assumption of detailed balance at the macroscopic scale, as ergodic conditions are assumed to hold only at small length scales. The derived REV condition is a constraint that fluctuations must obey to derive phenomenological equations that provide a valid representation for the entropy production in the system. Experimental data are then used to assess the role of capillary fluctuations in the context of the theory.

### A. Thermodynamic description

Thermodynamics provides a mechanism to link flow processes with conservation of energy. For non-equilibrium systems, this is accomplished by deriving an expression for the entropy production, which can then be exploited to derive phenomenological equations for particular processes.^{10} At the practical level, thermodynamics is a multi-scale book-keeping system to keep track of how energy is distributed within a system. The thermodynamics of heterogeneous systems are typically treated by dividing the system into sub-regions based on the system composition.^{12,45,46} In addition to the entropy *S*, we consider the extensive variables to be

*V*—fluid volume for $i\u2208{w,n}$,_{i}*A*—surface area for fluid_{n}*n,**A*—surface area for the solid material,_{s}*h*—film thickness along the solid material, and_{s}*N*—number of molecules of type_{ik}*k*within any region $i\u2208{w,n,s}$.

For simplicity energy associated with the solid is omitted, such as effects associated with deformation. Considering solid that does not deform, the surface area bounding the wetting fluid *A _{w}* is omitted based on the fact that only two of ${An,Aw,As}$ are independent. The internal energy is given as a function of the independent extensive quantities

^{71}

The Euler equation for the internal energy of the system is then written as

where the intensive measures are

*p*—fluid pressure for $i\u2208{w,n}$,_{i}*γ*—surface energy for fluid–fluid interface,_{wn}*γ*—surface energy for the solid material,_{s}Π

_{s}—disjoining pressure along the solid material, and*μ*—chemical potential for molecule_{ik}*k*within any region $i\u2208{w,n,s}$.

The indices indicate the entity where the property is defined. For example, the fluid pressures are defined over the fluid phases $i\u2208{w,n}$. The surface tension is defined over an interface region, which may be the interface between fluids, *wn*, or the fluid–solid interface *s*. In a strongly wetted system, the solid surface is coated everywhere by a thin film of fluid *w*. This means that fluid *n* interacts with *w* everywhere along its boundary, described by the surface energy *γ _{wn}*. If the system is not strongly wetted by

*w*, the formulation can still be used based on the fact that any excess energy due to interactions between the fluids and the solid will be accounted for by local variations in

*γ*along the solid. The complexity of the fluid–solid interactions is therefore embedded within the associated energy term, $\gamma sAs$. In our approach, we do not treat the solid as being in contact with a particular fluid. Instead, we consider the situation where the associated interfacial energy

_{s}*γ*and film thickness

_{s}*h*can vary along the solid. The formulation presupposes that films coat the solid surface everywhere, with the understanding that if the film thickness

_{s}*h*= 0 is locally on portions of the solid surface, this will be equivalent to treating films as covering only part of the solid surface. The disjoining pressure is defined as

_{s}The energy due to fluid–solid interactions is therefore lumped into two terms, one associated with the interfacial energy and a second due to the disjoining pressure and film contributions. Variations in *γ _{s}* and

*h*along the surface of the solid occur due to the interaction between the fluid and solid material, including the effects of roughness and chemical heterogeneity.

_{s}The terms that account for capillary fluctuations are neglected in conventional non-equilibrium thermodynamics; these appear only when time-averages are introduced. When considering flow processes that occur over long periods of time, it is easy to overlook fast dynamics. Haines jumps are one such category of event, with dynamics that typically occur over a duration of milliseconds to seconds.^{2,47,48} Topological changes such as droplet coalescence also dissipate energy, occurring even more rapidly and contributing to pressure fluctuations. While both classes of event are common for flows in porous media, the associated energy dissipation has not been properly included in thermodynamic models. In Secs. II B–F, we demonstrate that pressure fluctuations represent an essential contribution to the energy dynamics for two-fluid flows in porous media.

### B. Thermodynamic averages in space and time

The use of time averaging to study the non-equilibrium behavior for two-fluid flows in porous media was first treated by Aryana and Kovscek.^{49} While these authors demonstrated that the conventional Buckley–Leverett model could be derived based on this approach (as well as extended models), aspects pertaining to the non-equlibrium thermodynamics were not considered. Here, we develop this theory explicitly. Averages are constructed with the time-and-space averaging operator, defined as

where Λ is a region of time with duration *λ* and Ω is a spatial averaging region.^{50} $V$ is the volume associated with thermodynamic measurements which can be estimated from the length scale given in Eq. (2). For example, if the fluid pressure is measured with a transducer, $V$ can be thought of as the volume of fluid that is in local equilibrium with the measured value as a consequence of local mixing due to the motion of molecules in the immediate vicinity. The total volume of the spatial averaging region Ω is $V>V$. Dividing Eq. (7) by $V$ leads to a representation of the internal energy on a per-unit-volume basis. Since a point in space may be within one fluid or the other, we define the indicator function

where $i\u2208{w,n}$ to denote the region of space occupied by each fluid, Ω_{i}. Integrating $\varphi i$ over any region will therefore return the volume of the associated fluid within that region.

Statistical definitions for thermodynamic quantities rely on the ergodic hypothesis such that ensemble averages are the same as averages in space or time. The basis for this assumption is that the energy micro-states are locally well-mixed, i.e., that the molecules in the system interact and exchange energy on a timescale that is sufficiently fast relative to the size of the system considered. Within this context, Eq. (7) holds in the limit of infinite time and non-equilibrium expressions must be developed to study time-dependent processes. In the context of ergodic theory, the averaging operator given by Eq. (9) explicitly mixes information in space and time to generate averages. Non-ergodic effects can be accounted for in the averaged system provided that Λ and Ω are chosen appropriately. Averages of the extensive quantities are simple averages so that these quantities retain their basic meaning,

The intensive measures are then defined to maintain scale-consistency with Eq. (7),

The internal energy in the averaged thermodynamic model will then be given by

In a non-equilibrium system, deviations of the intensive quantities contribute to the energy dynamics.^{50} The deviation terms are defined as the difference between the microscopic value and the averaged value,

The time derivative for the internal energy includes contributions from the time derivatives of the averages of the extensive measures as well as fluctuation terms that arise due to the microscopic deviations in the intensive quantities,

where a sum is implied for all chemical potential terms. The multi-scale fluctuation terms account for transient effects and heterogeneity associated with the thermodynamic behavior:

$\u27e8Vi\u2202pi\u2032\u2202t\u27e9$—fluid pressure fluctuations, e.g., due to the transient effect of changing capillary forces on the fluid pressure during flow;

$\u27e8An\u2202\gamma wn\u2032\u2202t\u27e9$—fluctuations to the fluid–fluid interfacial tension, e.g., due to different surface tension along the fluid–fluid interface, which may arise if surfactants are present;

$\u27e8As\u2202\gamma s\u2032\u2202t\u27e9$—fluctuations to the fluid–solid surface energy, e.g., which occur due to movement of the contact line and complex fluid–solid interactions that occur due to surface heterogeneity;

$\u27e8Ashs\u2202\Pi s\u2032\u2202t\u27e9$—fluctuations to the disjoining pressure of films on the solid surface, e.g., which can occur due to localized film swelling; and

$\u27e8Nik\u2202\mu ik\u2032\u2202t\u27e9$—fluctuations to chemical potential, which arise due to mass diffusion.

There is a direct link between the pressure fluctuation defined in Eq. (15) and rugged energy landscapes depicted in Fig. 2. Pressure fluctuation terms contribute to the non-equilibrium description based on pore-scale events that transpire within the averaging window defined by Λ and Ω. The reversible work associated with an *ison* is captured based on the product of the fluid pressure difference and rate of volume change, as in standard non-equilibrium thermodynamics. During a strictly quasi-static displacement, the rate of change for the microscopic pressure *p _{i}* should closely correspond to the average fluid pressure $p\xafi$. On the other hand, if we consider a

*rheon*to be a rapid pressure change at constant volume, then

*p*must deviate significantly from $p\xafi$ over the duration of the event. Therefore, if there are

_{i}*rheons*occurring within Ω during the time interval Λ, then the associated fluctuation terms are expected to have a significant effect. This will be treated in more detail in sections to follow.

### C. Rate of work and entropy inequality

The change to the internal energy of the system is determined based on the rate of work done on the system and the heat added,

The instantaneous rate of work is given by the sum of Eqs. (3) and (5), the latter being already built into the thermodynamics. For the data considered in this work, which correspond to immiscible two-fluid flows in strongly water–wet systems, several assumptions can be made to simplify the form of Eq. (15):

- the system temperature is uniform and no heat is added to the system, which is justified by the large heat capacity for the system,$T=T\xaf\u2192\u27e8S\u2202T\u2032\u2202t\u27e9=0\u2009,\u2003\u2202Q\u2202t=0\u2009;$
- the fluid–fluid interfacial tension is constant everywhere,$\gamma wn=\gamma \xafwn\u2192\u27e8An\u2202\gamma wn\u2032\u2202t\u27e9=0\u2009;$
- the solid does not change shape, meaning that$\u2202A\xafs\u2202t=0\u2009,\u2003\u2202V\xafw\u2202t=\u2212\u2202V\xafn\u2202t\u2009;$
- compositional effects can be ignored,$\mu \xafik\u2202N\xafik\u2202t=0\u2009,\u2003\u27e8Nik\u2202\mu ik\u2032\u2202t\u27e9=0.$

These assumptions will typically hold for experiments designed to study immiscible fluid displacement in porous media. In surfactant systems, the fluid–fluid interfacial tension may vary within the system. Assumption 2 would not hold in such a system. In liquid–gas systems, effects due to compressibility may also contribute, which may lead to additional contributions from the chemical potential, i.e., the relaxation of Assumption 4. Subject to these assumptions and using Eq. (16) to replace the change to the internal energy with the rate of external work done on the system, Eq. (15) can be simplified and rearranged to obtain an entropy inequality,

where *L* is the distance between the sample inlet and outlet and the time-average for the fluid volume fraction is

In the context of Fig. 2, specific terms appear that correspond to the reversible pressure–volume work associated with an *ison* and the spontaneous *rheons*. Note that the average pressure will behave as a smooth function of time, meaning that the “rugged” parts of the energy landscape can be fully transferred to the fluctuation terms if Λ is large enough. Three fluctuation terms are retained, the first associated with capillary fluctuations in the fluid pressures. The second is due to wetting fluctuations that are attributed to the movement of the contact line during fluid displacement, which appears due to the time derivative for the deviation of the solid interaction energy $\gamma s\u2032$. Intuitively, the fluid–solid interaction energy changes when the contact line moves. Finally, the effects of film swelling contribute based on the fluctuation of the disjoining pressure, which will change based on changes to the film thickness. This contribution is important primarily when the wetting fluid is poorly connected, in which case films that coat the solid surface provide a primary flow pathway for that fluid. To obtain a more useful form of the entropy inequality, a flux–force form must be obtained.^{9} In Sec. II F, we demonstrate how this can be accomplished for systems where the fluctuation terms play a significant role.

### D. Fluctuations and solid wetting energy

The model considered here is developed for strongly water–wet systems where the solid is coated everywhere by a thin water film. This means that the boundary of the non-wetting fluid touches water everywhere, and the surface energy is proportional to the total boundary surface area *A _{n}*. This simplifies the geometric analysis presented in Sec. II E. For the situations considered here, where the surface wetting properties are homogeneous and do not change with time, the surface wetting fluctuation is zero,

In many situations, the surface wetting properties will be heterogeneous, with the local wetting energy varying based on chemical properties and roughness of the solid. In these more general conditions, all of the solid surface energy contributions are embedded into a single term. This term can also account for situations where the solid is not strongly wetted by one fluid. In this case, only part of the solid is in contact with water, with associated surface area *A _{ws}*. The rest of the solid contacts the non-wetting fluid, with $Ans=As+Aws$. In this situation, the fluctuation term can be subdivided into distinct regions of water and non-wetting fluid contact of the grain surface,

If the solid surface energy is constant over each region with $\gamma s=\gamma \xafws$ everywhere on *A _{ws}* and $\gamma s=\gamma \xafns$ everywhere on

*A*, then

_{ns}In this case, fluctuations of the solid surface energy are determined by movement of the contact line, which changes the surface areas. Since $\gamma \xafns,\u2009\gamma \xafws$, and $A\xafs$ are constant with respect to time,

Since the total solid surface area is constant

Using this result with Eqs. (21) and (22) and doing some basic algebra show that the fluctuation term is

This is the conventional expression for change in fluid–solid surface energy. Note that while *A _{ns}* and

*A*are extensive quantities, only their sum,

_{ws}*A*was included in Eq. (6). Since

_{s}*A*and

_{ns}*A*are identified explicitly from heterogeneity in an intensive property,

_{ws}*γ*, all associated energy will be embedded into the fluctuation term for that intensive property. Note that Eq. (24) is a simplification of the more general case, where

_{s}*γ*may vary in complex ways along the grain surface, or change with time due to transient chemical effects. We now consider particular geometric considerations for water–wet systems.

_{s}### E. Geometric constraints

As the fluid is injected into a porous material, the configuration of fluids within the micro-structure changes, with interface rearrangement occurring alongside displacement. Changes to the fluid volume fraction will typically also lead to changes in interfacial area. However, since the interfacial area can change without any change to the volume fraction, the chain rule cannot be applied to simplify Eq. (17),

The rate of change in surface area is neither independent from the rate of change in the volume fraction nor entirely dependent on it. To obtain a flux–force form of the entropy inequality, the dependent component must be distinguished from the independent component. It has been previously showed that a global geometric relationship can be constructed based on a non-dimensional relationship developed based on invariant measures of integral geometry.^{24} Here we apply this relationship to the time-and-space averaged geometric measures,

The non-dimensional measures are defined as

where $A\xafn$ is the time average for the surface area of the non-wetting fluid boundary, $H\xafn$ is the time average of integral mean curvature, and $\chi \xafn$ is the time average for the Euler characteristic. Even though the instantaneously measured $\chi n$ is a discontinuous function of time, $\chi \xafn$ will be continuous due to its construction as a time average. While the Gibbs dividing surface can be considered as introducing non-smooth aspects into the system description based on the discrete representation of fluid regions according to set theory, the time-and-space average removes any non-smooth behaviors by incorporating the past and future state of the system within the instantaneous representation. To link the time derivatives for the surface area and volume fraction, we note that a differential form is implied by Eq. (25),

Since $W\xafn$ and $X\xafn$ include $A\xafn$, this can be arranged to remove the rate of change in surface area from Eq. (17). The link between thermodynamics and the geometric evolution can be understood by expanding the derivatives using the definitions in Eq. (26),

Topological contributions arise based on the time derivative of $X\xafn$,

Inserting the previous two expressions into Eq. (27) and rearranging terms lead to an expression that links the time derivatives,

where three geometric functions have been determined to depend on the geometric state of the system,

The average curvatures for the surface can be obtained by dividing the integral mean curvature and Euler characteristic by the total surface area,

This means that

where the average mean curvature $J\xafn$ and average Gaussian curvature $K\xafn$ are now cast as intensive properties of the surface. Inserting into Eq. (30) and rearranging terms lead to

and expanding terms, we find that

and

which means that

Solving to eliminate the time derivative for the surface area as it appears in the thermodynamic expressions,

### F. Flux–force form of entropy inequality

Classical non-equilibrium theory relies on a flux–force form of the entropy inequality to derive phenomenological equations.^{10} Using the fact that the molecular system is time-reversible, detailed balance is assumed.^{9} However, in systems with cooperative effects, such as interface rearrangements, detailed balance will not hold within sub-regions of the system occupied by individual fluids. The fluctuation terms are of central importance to this effect, since they capture the net contribution due to cooperative effects over Λ and Ω. As long as the fluctuations do not undermine the independence of terms in the flux–force entropy inequality, phenomenological equations may be derived. Using the results of Secs. II A–E, we now obtain such a result. We consider situations where changes to the film thickness are caused by changes to the fluid volume fraction,

Making use of Eqs. (41) and (42), the entropy inequality in Eq. (17) can be put into a flux–force form

Provided that Ω and Λ are sufficiently large, the averaged quantities will be smooth functions of time, since these are time-averaged quantities. The fast dynamics will be entirely embedded within the associated fluctuation terms. Equation (43) is not yet in a flux–force form due to the contribution from fluctuations and the geometric evolution terms. Physically, this aligns with transient redistribution of energy within the system. For example, there may be reversible energy transfer from the pressure fluctuation terms to the surface energy. Note that a flux–force form can be obtained if

This expression defines a condition for the representative elementary volume (REV) for two-fluid flow through porous media. If terms on the left-hand side of Eq. (44) do not sum to zero, then they constitute a net contribution from fluctuations to the energy dynamics. If this happens, terms on the first two lines of Eq. (43) will not be independent due to couplings caused by the internal energy dynamics. Assessment of terms in Eq. (44) is therefore fundamental to the validity of macroscopic equations for two fluid flow.

To evaluate Eq. (44) in experiments, pressure fluctuations need to be evaluated for both wetting and non-wetting fluids. The geometric state variables are needed to calculate the time derivatives for curvature terms and the slope of the geometric state function. We note that the fluctuations do not necessarily need to be Gaussian to obtain a flux–force entropy inequality. As long as any asymmetry in the fluctuations is balanced by changes to the surface energy based on the evolution of interface configuration, Eq. (44) can be satisfied. Fluctuations should be considered from two perspectives: (1) steady-state flow at constant fluid volume fraction; and (2) unsteady flow where the fluid volume fraction is changing. In a steady-state flow, the curvature must be independent of time, since any net change to the geometric configuration of fluids is inconsistent with the system being at steady-state. Therefore, the associated conditions are

There must be no net rate of energy change from fluctuations at steady-state, since net changes to the system are inconsistent with the basic notion of steady state. For unsteady displacement, capillary fluctuations may be balanced by changes to the surface energy. This reflects the fact that the fluid configuration can change at constant fluid volume fraction based on changes to the boundary curvature. Since the boundary curvature can evolve independently of the volume, a separate term must appear to account for these effects.

Since the terms in Eq. (44) do not multiply the change in fluid volume fraction in Eq. (43), it is essential that their average value be zero. Otherwise, the flux–force form of the entropy inequality cannot be exploited to derive phenomenological equations. Any net contribution from the fluctuation terms during steady-state flow invalidates any corresponding measurement of the relative permeability. On the other hand, if Eq. (44) is satisfied, then the effective permeability coefficients can be considered to fully account for the dissipation, including contributions from pore-scale events on the basis that the average rate of events is constant based on the averaging domain, Λ and Ω. Results obtained using pore-network models for a steady-state fractional flow suggest that the associated fluctuations will be symmetric.^{7} In this case, a flux–force form is possible for Eq. (43), with the thermodynamic forces identified as

The resulting flux–force form of the entropy inequality can naturally be applied to derive the standard constitutive relationship for relative permeability, given in Eq. (4). This conventional relationship assumes that the relationship between the flow rate $q\xafi$ and the applied forces is linear, which is known to break down as the capillary number increases. This result is also based on the assumption that the forces and fluxes are independent, which has been called into question based on the results of homogenization theory.^{13,51,68} Previous efforts to derive Eq. (4) from first principles have generally concluded that the flow rates $q\xafw$ and $q\xafn$ are not independent due to momentum exchanges between the two fluids. For example, if an external potential gradient is applied only to one fluid, it will pull some amount of the other fluid along with it, which is not accounted for based on Eq. (4). Cross-coupling forms have been proposed to account for this effect.^{51,69} However, since cross coupling forms introduce additional phenomenological parameters, it is considerably more complicated to measure experimentally and therefore difficult to apply in real systems. From the more practical standpoint, one can consider relative permeability as being a proxy for the actual dissipation rate, which depends non-linearly on the fluid geometry, the mobility ratio capillary number, and other factors.^{53–55} Since entropy is a scalar, we can always embed this rate within a scalar parameter. More recent efforts align with this more practical conceptual picture and can be considered as characterizing the non-linear dependence of the actual dissipation rate on the applied forces.^{52,70} In this context, the fluctuation constraint in Eq. (46) is a criterion that must be satisfied to accurately measure the steady-state dissipation, since this ensures that the energy dynamics are stationary.^{68,70,71}

The rate of change for the fluid volume fraction has been previously used to derive expressions for capillary pressure dynamics.^{11,44} We obtain an alternate form that includes the effect of both films and capillary fluctuations,

where *τ _{w}* is a phenomenological parameter. Equation (48) is an expression for the capillary pressure dynamics, which differs from previous forms. First, the effect of film swelling is included, which is important for the snap-off mechanism. The capillary term includes the contribution over the entire fluid boundary, including the part of the boundary where the shape is strongly influenced by the local solid micro-structure. Due to the effect of films, the force-balance in these regions differs from the meniscus. In the averaged form, the equilibrium requirement is that all force contributions must cancel, including the effect of the fluid pressures, meniscus curvature, and film contributions. Averaging in both time and space ensures smoothness of the pressure signal with respect to time provided that the system size and time interval for averaging are sufficiently large. Second, the capillary fluctuation terms are accounted for explicitly based on the REV constraint stated in Eq. (44). This REV constraint is the main assumption used to derive Eq. (48). It is possible that net contribution due to capillary fluctuations will vanish in REV-scale systems due to the cancelation of terms. From the perspective of pore-scale events, Ω and Λ must be large enough to capture the distribution of events within a particular material. This will depend on both the material and the flow rate. However, the fluctuation term is not guaranteed to vanish because the pressure fluctuation signal may be non-stationary during unsteady displacement. If this is the case, then Eq. (43) does not reduce to a flux–force form and additional treatment would be required to derive macroscopic phenomenological equations. It is therefore important to treat the fluctuation terms cautiously. In Sec. III, we consider several practical examples to illustrate the behavior for capillary fluctuations and their contribution to the energy dynamics in multiphase systems.

## III. RESULTS

### A. Flow through porous media

The flow of immiscible fluids in porous media is strongly influenced by geometric and topological effects.^{56–64} Capillary forces typically dominate such flows, which are directly impacted by changes to fluid connectivity. We consider fluid displacement within a porous medium as a means to explore how geometric changes influence physical behavior. A periodic sphere packing was used as the flow domain. Using a collective rearrangement algorithm, eight equally sized spheres with diameter *D* = 0.52 mm were packed into a fully periodic domain $1\xd71\xd71$ mm^{3} in size. The system was discretized to obtain a regular lattice with 256^{3} voxels. Simulations of fluid displacement were performed by injecting fluid into the system so that the effect of topological changes could be considered in detail. Initially the system was saturated with water, and fluid displacement was then simulated based on a multi-relaxation time (MRT) implementation of the color lattice Boltzmann model.^{66}

To simulate primary drainage, non-wetting fluid was injected into the system at a rate of *Q* = 50 voxels per time step, corresponding to a capillary number of $4.7\xd710\u22124$. A total of 500 000 timesteps were performed to allow fluid to invade the pores within the system, leading to the sequence of structures shown in Figs. 4(a)–4(c). Fully water–wet conditions were prescribed, with the viscosity and density ratio both equal to one, representing a generic displacement involving dense fluids. The simulation conditions were constructed to demonstrate topological changes in a relatively small system, including both loop formation and snap-off. The first pores are invaded at *t* = 0.19 pore volumes, which is just after the system overcomes the entry pressure. The fluid does not yet form any loops, since a pore must be invaded from two different directions before a loop will be created. Haines jumps and loop creation can be considered as distinct pore-scale events that may coincide in certain situations. At time *t* = 0.34 pore volumes, the first loop is formed. As primary drainage progresses additional loops form, causing the Euler characteristic tends to decrease as non-wetting fluid is injected.

At the end of primary drainage, the flow direction was reversed to induce imbibition, with water injected at a rate of *Q* = 20 voxels per time step. Sequences from the displacement are shown in Figs. 4(d)–4(f). As water imbibes into the smaller pores in the system, there is a well-known tendency to snap off non-wetting fluid by destroying connections to larger adjacent pores where imbibition has not yet occurred. This is known as the Roof snap-off mechanism.^{39} Snap-off breaks loops within the system, causing the Euler characteristic to increase. The order that loops snap off during imbibition is not the reverse of the order that the loops form during drainage. This is because the drainage process is dominated by the size of the throats whereas imbibition is governed by the size of the pores. Eventually, the sequence of snap-off events can cause non-wetting fluid to become fully disconnected from the main region, leading to trapping as seen at time *t* = 1.86 pore volumes.

The presence of trapped fluid distinguishes secondary drainage from primary drainage. The sequence of fluid configurations shown in Figs. 4(g)–4(i) depicts the sequence of coalescence events that reconnect trapped fluid. Due to the presence of trapped fluid, the order that loops are formed within the porespace does not match the order observed during primary drainage. Primary and secondary drainage each show a rapid spike in the pressure difference due to the pore entry pressure. The entry sequence along primary drainage matches the configurations from Figs. 4(a)–4(c) and corresponds to events labeled as A and B in Fig. 5(a). The high pore entry pressure is consistent with expected behavior that the largest energy barriers in the capillary dominated system are due to pore entry.^{2} The pore entry sequence along the secondary drainage is shown in Figs. 4(g)–4(i) and labeled as D and E in Fig. 5(a). On secondary drainage, the entry pressure is slightly reduced due to the presence of trapped fluid. This shows that the first coalescence event depicted in Figs. 4(g)–4(i) corresponds exactly to the maximum fluid pressure difference observed during pore entry. This is distinct from the initial pore entry observed during primary drainage; upon secondary drainage a fluid coalescence event occurs instead of a Haines jump. The resulting pressure fluctuations differ as a consequence.

Like coalescence, snap-off events are associated with clear fluctuations in the bulk fluid pressures. The snap-off sequence shown in Figs. 4(d)–4(f) is labeled as C in Fig. 5(a). As loops of fluid are destroyed, corresponding fluctuations in the fluid pressure difference are observed. The largest disruption is associated with the snap-off event. The capillary pressure for the connected and disconnected regions of non-wetting fluid is shown in Fig. 5(b). Consistent with the observation of other authors, fluid is trapped at a higher capillary pressure than what is measured in the connected fluid phases at the instant of snap-off.^{39} As soon as snap-off occurs, the time derivative of the pressure difference between the two connected fluid phases immediately jumps from strongly negative to strongly positive. This occurs as each fluid region relaxes toward its preferred equilibrium state, which is possible because the two fluid regions are no longer mixing after the snap-off event. The final pressure for the trapped region is determined based on geometry, with the fluid assuming a minimum energy configuration within the pore structure where it is trapped.

The energy associated with pressure fluctuations that arise due to pore-scale events is shown in Fig. 6. At the pore level, the energy due to fluctuations is comparable to the rate of pressure volume work and the rate of change in surface energy. Defining the exchanges between these terms in a non-equilibrium setting is non-trivial, since it is determined in part by the geometry of the problem. It is clear from these results that the fluctuation terms contribute a significant amount to the overall energy dynamics at the scale of an individual pore. We now extend this analysis to larger scale systems where large numbers of pore-scale events occur, allowing their frequency and distribution to be analyzed in detail.

### B. Experimental results

Experimental data were collected for a Berea sandstone.^{2} The porosity of the sample was measured using the Archimedes method based on a physical measurement for a one inch core plug. This particular Berea block had an extremely homogeneous structure and from tens of samples the relative standard deviation was just 0.2% in porosity units.^{67} The sample used for imaging had a diameter of 4 mm and length 10 mm. The sample porosity was 0.199 and the absolute permeability was 0.7 *μ*m^{2}. Oil was injected into the sample at a rate of 0.35 *μ*l/min. Pressure transducer measurements were collected for the oil at an interval of $\Delta t=0.32$ s. Several assumptions were necessary to support the analysis of the experimental data. Since pressure measurements were available for the oil phase only, the water was assumed to be in equilibrium with the external pressure (1 atm) at the outlet of the sample. Since the transducer measurements for the oil pressure are available only at the inlet boundary, this was assumed to adequately represent the spatially averaged pressure within the sample. While it is nearly certain that various intermittent pressure gradients and spatial fluctuations occurred within the sample, the response of the pressure waves (based on the speed of sound) is fast compared to the timescale for flow (based on the rate of change in saturation). Time derivatives were evaluated by fitting a spline function to the experimental time series data and evaluating the first derivative. Time averages were then computed based on a simple moving average (SMA). Average pressure values are shown in Fig. 7. For Λ = 128 s, the pressure signal behaves as a continuous *ison*, incorporating the reversible components of the original signal. Given a sufficiently long time averaging interval, the effect of *rheons* is thereby embedded entirely within the fluctuation terms.

Figure 7 shows the experimental pressure signal and time-averaged curves computed based two time intervals, Λ = 16 s and Λ = 128 s. The associated fluctuation distributions are shown in Fig. 8. Pressure fluctuations continue to be evident due to the larger events with Λ = 16 s. These fluctuations are mostly removed for Λ = 128 s, although the effect of larger events remains visually evident. Time average effectively removes the noise due to the fast pressure fluctuations from the pressure signal. The distribution for the pressure fluctuation is shown in Fig. 8. Since the pressure fluctuation is defined based on the difference between the instantaneous and average pressure values, the amplitude is sensitive to the averaging interval. Figure 8(a) shows the distribution for Λ = 16 s. The distribution for Λ = 128 s is shown in Fig. 8(b). According to the central limit theorem, the distribution should be well-approximated by the Gaussian distribution provided that Λ is large enough to include a sufficiently large number of pore-scale events. Standard assumptions also require that the average pressure $p\xafi$ be stationary; if the average pressure value is drifting in time, then skew may be introduced in the pressure fluctuation distribution. This is an effective condition for the REV of the system. As the system grows larger, the rate of change to the fluid volume fraction $\u2202\varphi \xafw/\u2202t$ will be slower, since it takes a longer time for flow to fill the larger system. However, the timescale for pore-scale events will be the same. The average pressure signal will tend to be more stationary for larger systems considered over longer time intervals. For a fixed system size, a larger Λ will smooth the pressure signal and reduce the amplitude of the pressure fluctuations, as shown in Fig. 9. However, a larger value of Λ will mean the distribution is less stationary, since the average pressure can drift more over a longer time interval. For example, it is clear from Fig. 7 that the average pressure difference is changing at the timescale for Λ = 128.

A reasonable criterion for the REV is that the magnitude for the pressure fluctuations (based on the variance for the pressure fluctuation distribution as depicted in Fig. 8) must be small as compared to the rate of pressure–volume work. If this condition is not met, then *rheons* will exert an asymmetric influence on the pressure signal and Eq. (44) may not hold. In this situation, the average energy landscape is rugged in nature and the fluctuations must be considered explicitly. In a sufficiently large system, considered over a sufficiently large time interval, Eq. (44) should hold and an averaged non-equilibrium description should be possible. Figure 9 shows the amplitude of pressure fluctuations on the energy scale, which includes the rate of pressure–volume work. The instantaneous values for the pressure fluctuation can be very large, due to the fact that pore-scale events are very fast. The overall energy contribution is significant and cannot be ignored unless REV requirements are met. An intriguing possibility is that various terms in Eq. (44) may counter-balance asymmetry in the pressure fluctuations. For example, pressure fluctuation energy can be redistributed to the surface energy, or to the other fluid. Part of the energy from *rheons* can be reversibly transferred to other internal energy modes.

An assessment of the temporal REV is shown in Fig. 10. Assessment of the REV based on the pressure fluctuations must consider two essential components: (1) the mean for the pressure fluctuation distribution, since a non-zero mean implies there is a net contribution to the internal energy; and (2) the variance associated with the amplitude of the pressure fluctuation, which is related to the number of pore-scale events within the time interval Λ. For the particular experimental data considered here, the energy associated with the mean pressure fluctuation is $\u223c5$% of the rate of pressure–volume work. This can likely only be further reduced by increasing the physical size of the sample. The standard deviation for the pressure fluctuation decreases as Λ increases, decreasing to $\u223c10$% of the rate of pressure–volume work for Λ = 512 s. When considered over short time averaging intervals, the fluctuations are several times the rate of pressure volume work. Based on the metrics defined in Fig. 10, the REV accuracy for the experimental system considered is likely to be in the range of ±5%–10% for the capillary pressure.

It is noted that the requirement given in Eq. (44) can still be satisfied if the net contribution from the pressure fluctuation balances against the change in surface energy and the contribution due to film swelling. However, this cannot be assessed for the data set considered here, since the data needed to evaluate the curvature and film contributions were not collected. Additional study of this relationship is important to establish the validity of macroscopic flow equations. Our results suggest that during unsteady flow, the pressure fluctuations within an individual fluid will make a net contribution to the energy dynamics that cannot be ignored. To meet an REV criterion, these fluctuations must either cancel with other terms in Eq. (44), or the net effect of fluctuations must vanish for sufficiently large systems.

It is important to recognize that pore-scale events will always dissipate a significant amount of energy during two-fluid flows in porous media. Averaging does not change this. What averaging does accomplish is to treat fluctuations such that their energy dissipation is accounted for based on Eqs. (4) and (48), i.e., in situations where Eq. (44) is satisfied. The constitutive models for relative permeability and dynamic capillary pressure are not valid if these constraints are not met. However, when the system size is large and longtime intervals are considered, the net contribution of the capillary fluctuation terms will trend toward zero. Their role on dissipation within the system will then be incorporated into the phenomenological coefficients *k _{rn}*,

*k*, and

_{rw}*τ*.

_{w}## IV. SUMMARY

We derive non-equilibrium thermodynamics for immiscible fluid flow in porous media using the method of time-and-space averaging. Averaging in both time and space reveals macroscopic thermodynamic fluctuation terms, which are athermal in nature and account for internal energy transfers due to multi-scale effects such as capillary pressure dynamics. An entropy production equation is derived that includes the relevant fluctuation terms, which is converted into flux–force form to support the derivation of phenomenological equations. The standard equations for effective permeability can be derived, which can be generalized to cross coupling forms in situations where significant momentum exchanges occur between fluids. Key results of the theory are (1) a time-and-space REV constraint on fluctuations observed during steady-state displacement, which is needed to derive the flux–force form of the entropy inequality; and (2) a dynamic capillary pressure expression that accounts for the effects of films and capillary fluctuations. Using simulation and experimental data, we show that capillary fluctuations represent a significant part of the non-equilibrium energy dynamics for immiscible fluid flows. Based on experimental data, we assess the distribution for capillary fluctuations for time intervals ranging from Λ = 8 s to Λ = 512 s. The net contribution for fluctuations during unsteady displacement varies from 5% to 10% of external pressure–volume work. Our results indicate that capillary fluctuations are an essential consideration to determine REV for immiscible displacement. Further work is needed to determine whether net energy contributions due to capillary fluctuations and surface energy will cancel based on the REV constraint derived in this paper.

## ACKNOWLEDGMENTS

An award of computer time was provided by the Department of Energy Summit Early Science program. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract No. DE-AC05-00OR22725. *μ*CT was performed on the TOMCAT beamline at the Swiss Light Source, Paul Scherrer Institut, Villigen, Switzerland. We are grateful to G. Mikuljan at Swiss Light Source, whose outstanding efforts have made these experiments possible'.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## References

_{2}sequestration