The Relation Between Dissipation and Memory in Two-Fluid Displacements in Disordered Media

) framework to sort out and evaluate the energy dissipated in the passage between metastable equilibria. This framework provides a link between microscale quantities (microscopic quenched disorder, metastable equilibrium configurations) and the upscaled macroscopic trajectories and energy dissipation. The procedure is generic; here we apply it to a model of fluid-fluid displacements in an open fracture Abstract We show that the return-point memory of cyclic macroscopic trajectories enables the derivation of a thermodynamic framework for quasistatically driven dissipative systems with multiple metastable states. We use this framework to sort out and quantify the energy dissipated in quasistatic fluid-fluid displacements in disordered media. Numerical computations of imbibition–drainage cycles in a quasi-2D medium with gap thickness modulations (imperfect Hele-Shaw cell) show that energy dissipation in quasistatic displacements is due to abrupt changes in the fluid-fluid configuration between consecutive metastable states (Haines jumps), and its dependence on microstructure and gravity. The relative importance of viscous dissipation is deduced from comparison with quasistatic experiments.

HOLTZMAN ET AL. 10.1029/2023GL104073 2 of 8 (an imperfect Hele-Shaw cell) where the microscopic physical mechanisms of surface tension and capillarity lead to macroscopic pressure-saturation trajectories exhibiting the complex behavior of hysteresis and RPM.

Generic Formulation for Fluid-Fluid Displacements
Energy conservation ensures that between consecutive metastable equilibria the mechanical work invested in driving the fluids is partially stored as internal energy of the multiphase system configuration, and partially dissipated. Thus, for an infinitesimal change in wetting-phase saturation (corresponding to a volume change dV w ), the energy balance applies, where we have adopted the convention that đW > 0 is the external mechanical work performed on the system to drive the displacing fluid (đW < 0 if extracted; the sign depends on the external force and direction of advancement), dU > 0 is the increase in internal energy of the two-fluid configuration, and đΨ ≤ 0 is the amount of energy dissipated. The terms in Equation 1 refer to changes between two consecutive metastable equilibria, and đ denotes differentials of magnitudes that are not state functions (i.e., that depend on the path in the pressure vs wetting-phase saturation (PS) space). đΨ could be cast in terms of entropy production, but thermal fluctuations here are insufficient to bring the system over the large energy barriers between neighboring equilibrium configurations.
We consider first a generic scenario in which a fluid is injected at one side (inlet) of a disordered domain and displaces a second fluid. A specific example is the system shown in Figure 1. For a given disorder realization, the total energy of the system (its Hamiltonian  ) depends simultaneously on the interfacial configuration {ξ}, for example, the location and shape of the fluid-fluid interface (a microscopic feature), and the applied (macroscopic) pressure P on the fluid at the inlet; thus, every equilibrium state is defined by {ξ} and P. To obtain U and W, we split the Hamiltonian through the Legendre transformation  = − , where V w is the volume of the wetting phase in the domain (saturation times domain volume). In contrast with  , the internal energy U of the multiphase configuration depends only on the set of variables {ξ} and not on the sequence of driving pressures P. The crucial point comes now: if the RPM property holds, the configuration {ξ} is exactly recovered in a cyclic excursion of P. Thus, U = U({ξ}) also returns to its original value in a cyclic excursion, making it a true state function. We note that the internal energy of the compartment models developed by Cueto-Felgueroso and Juanes (2016) and Helland et al. (2021) is also a state function, thanks to the RPM property. In contrast, the amount of work depends on the path in the PS space, Figure 1. Schematic of the model and its experimental realization. Disorder is introduced by randomly placed "mesa defects"-sharp variations in gap thickness, b(x, y) = b 0 ± δb(x, y). Imbibition and drainage are driven by raising and lowering the wetting fluid reservoir, H. The cell is tilted to provide an effective gravity g e = g sinα. The interface configurations h(x) are recorded by a camera.
HOLTZMAN ET AL. (2) This expression relies on the incompressibility of the more wetting fluid, which guarantees that the amount of fluid displaced at the boundary of the domain coincides with the change of its volume in the medium. In this convention, the sign of đW depends on that of P (<0 in tension and >0 in compression) and dV w (>0 in imbibition and <0 in drainage).
Under quasistatic driving, a trajectory in the PS space is made of consecutive displacements of two types. In isons, the system remains trapped in a local energy minimum while P evolves exceedingly slowly and causes a smooth evolution of V w . In rheons, an abrupt change in V w occurs at constant P. Rheons take place when the limit of metastability of a local energy minimum is reached, the minimum disappears, and the system jumps to a new minimum, in what is called a Haines jump. In the idealized case of quasistatic driving, the timescales of these two kinds of displacements are infinitely separated, such that Haines jumps are effectively instantaneous in the time scale of P. This is the prototypical framework of spatially-extended athermal systems that undergo collective rearrangements (avalanches) under quasistatic driving (Jensen, 1998;Leschhorn et al., 1997;Pruessner, 2012). The dissipation between two consecutive equilibrium configurations (t − 1 and t) is obtained by integrating Equation 1, Ψ t−1→t is nonzero in rheons.
Also, since U is a state function, ∮dU = 0 in a closed PS cycle. Integrating Equation 3 along the cycle demonstrates that the total energy dissipated is the area encompassed within the cycle on the PS plane, Furthermore, since by definition Ψ cyc ≤ 0, there is only one sense allowed for contouring the cycle: for a given V w the drainage path occurs at lower P than the imbibition path-in agreement with experimental observations (Albers, 2014).

Model System: An Imperfect Hele-Shaw Cell
To gain further quantitative understanding of the microscopic mechanisms for dissipation, and enable rigorous comparison with experiments, we derive explicit expressions of W, U, and Ψ for quasistatic two-fluid displacements in an imperfect Hele-Shaw cell, in the framework of the model introduced by Holtzman et al. (2020).
The model considers the 2D projection of a heterogeneous Hele-Shaw cell, subjected to an effective gravity g e in order to prevent viscous fingering during drainage ( Figure 1). Gravity also allows investigating a wider range of pressure-saturation values than for a horizontal cell, as in the absence of heterogeneity there is no equilibrium position without gravity (Ayaz et al., 2020;Holtzman et al., 2020).
Disorder is provided by fluctuations in capillary pressure at the fluid-fluid interface, corresponding to localized modulations in the gap thickness, b(x, y) ("defects" forming constrictions and expansions). Fluid displacements are driven by an external pressure P = ρgH applied at the inlet, with ρ being the wetting fluid density, g the gravitational acceleration, and H the external head. The density of the second fluid is considered negligible (air).
Multiphase configurations {ξ} in the projected 2D model correspond to metastable equilibrium interface positions, {h(x)}. For a given disorder realization, they can be resolved from the pressure balance (Holtzman et al., 2020) The first term is the linear approximation (considering small interface deformations, | ℎ∕ | < 1) of the in-plane contribution to the Young-Laplace pressure jump across the interface, with γ the interfacial tension; the second is the hydrostatic pressure of the wetting fluid column at position x; the third is the external forcing; and the fourth p c (x, y) = 2γ cos θ/b(x, y) is the out-of-plane contribution to the Young-Laplace pressure jump at the front position h(x), with θ the apparent contact angle. Each disorder realization is defined by the microstructure b(x, y).
HOLTZMAN ET AL.
10.1029/2023GL104073 4 of 8 Considering mild local interface deformations we rule out the possibility of overhangs or complex processes such as snapoff.
The condition of mechanical equilibrium in Equation 5 can also be expressed as [ℎ( )] = − ∕ ℎ( ) = 0 , with the Hamiltonian given by This is equivalent to the Hamiltonian of a continuous Random Field Ising Model (RFIM) where p c (x, y) plays the role of the random field (Ganesan & Brenner, 1998;Grinstein & Ma, 1983). The terms in  as well as W, U, and Ψ in the forthcoming analysis are in units of energy per length, that is, in our 2D model they are normalized by the average gap thickness b 0 . Similarly, in 2D we use the area covered by the wetting fluid, , instead of its volume V w . Its relation to the wetting phase saturation * is given in Supporting Information S1. We note that for very large thickness variations, the integrals computing the energies should be corrected to take into account these variations.

Introducing Equation 6 into
 = − provides the internal energy (per unit thickness) of a given equilibrium configuration h(x), L is the width of the cell in the x direction. This expression of U accounts for the capillary energy of the in-plane and out-of-plane front deformations (first and third terms), and the gravitational potential energy of the wetting fluid. The work in Equation 2 for every elementary step dS w is đ The internal energy U is the energy of static equilibrium configurations, when the meniscus is at rest and the energy depends only on capillary and hydrostatic forces. According to Equation 3, the change in internal energy between two consecutive equilibrium configurations does not equate the work provided by the external driving; the difference is the energy dissipated, Ψ t−1→t , which in this framework is a loss of interfacial energy. This formulation does not account for viscous losses (viscosity of the fluids is not accounted for), in the expectation that they will be comparatively small in quasistatic displacements, as shown later.

Dissipation in Avalanches and the Role of System Properties
Dissipation depends on the interactions between the interface and the disordered medium. Surface tension introduces correlations along the interface, and the medium spatial heterogeneities turn these correlations into complex collective behavior (Holtzman et al., 2020). Spatial interactions among multiple defects of variable properties lead to a wide range of avalanche sizes (jumps in saturation) and dissipated energies.
To exemplify these properties of energy dissipation we compare the probability density function (PDF) of the energy dissipated within individual jumps (between adjacent equilibrium configurations t − 1 and t) normalized by the work input. We generate the disordered medium by randomly placing defects (Figure 1) of various strengths p c , drawn from different distributions (here, Gaussian with standard deviation and dichotomic; see Supporting Information S1). To establish the role of the disorder strength, we compare media of variable width of the defect strength distribution, (keeping the spatial distribution identical; see Supporting Information S1). Increasing , which sets the magnitude of strongest defects, increases the magnitude of the dissipation events, stretching the PDF toward larger values (see Figure 2a). Figure 2 also demonstrates that the energy dissipated in a single jump could greatly exceed the work invested in driving the system. This is because the energy released in an avalanche could have been stored as internal energy during many previous elementary steps. Most of the dissipation during a complete imbibition-drainage cycle therefore occurs in a handful of large events (Movies S1 and S2). This behavior is akin to the sudden release of energy in earthquakes (Sornette & Sornette, 1989) or granular avalanches (Denisov et al., 2016).
Another nonintuitive result is that the avalanche size (change in saturation) is not necessarily proportional to the amount of energy dissipated (Figures 2b and 2c). This non-proportionality is a general property of 5 of 8 quasistatically-driven disordered systems (Ortín & Goicoechea, 1998 , where 0 = 2 cos ∕ 0 and = 0 ( ∕ 0)∕(1 − ∕ 0) . The slope was found by fitting a straight line (in dashed gray) to the conditional average of the data (solid black line in insets of Figures 2b and 2c). The larger slope reflects the increased dissipation per avalanche size, due to the depinning from the stronger defects. We note the larger spread as well as deviation from a linear fit of smaller values, emphasized in log-log plots (Figures 2b and 2c insets).
Next, we examine how the system properties affect the total dissipation along a closed drainage-imbibition cycle, Ψ cyc . We find that Ψ cyc scales with the disorder strength ; for the Gaussian distributions studied here, Ψcyc ∼ ( ) with n ≈ 1.4 (Figure 3a). This nonlinear relationship is a result of interactions among the heterogeneities due to the lateral correlation caused by interfacial tension (Holtzman et al., 2020). Dissipation is also controlled by the effective gravity g e which sets the Bond number, namely the ratio of gravitational to capillary forces. Increasing g e reduces the capillary rise (Jurin's height), restricting the size of the avalanches | * , − * , −1 | and hence the dissipated energy, cf. Figure 3b and Movies S3 and S4.

Experiments Exposing Viscous Dissipation and the Limitation of the Quasi-Static Concept
The simulations above provide a quantitative analysis of the energy losses associated with the dissipation of interfacial energy (stored in the deformed interface configurations (Morrow, 1970)) in ideally quasistatic conditions.  The simulated PS trajectories show good agreement with an experiment in a disordered Hele-Shaw cell with similar material and geometrical properties, as shown in Figure 4. The experimental setup, procedure and analysis are detailed in Supporting Information S1. For the specific set of parameters in Figure 4 (dichotomic gap spacing b 0 = 0.46(1) mm and δb = 0.06(1) mm, disorder units of size 0.40(1) × 0.40(1) mm 2 covering 35% of the total area, and g e = 0.86(1) m/s 2 , see Supporting Information S1), the dissipated energy computed from the area enclosed within the PS cycle in the experiments is larger only by ∼18% than in the simulations. This implies that capillary losses account for most of the energy dissipated in quasistatic fluid displacements. Moreover, the primary imbibition curve and the early stages of subsequent drainage agree well, but the experimental data depart from the simulated drainage curve as the external pressure is lowered further. This also produces a small shift in the internal cycle. The lower external pressure required to drain a given wet area in the experiments, corresponding to a larger external work, is responsible for the larger dissipation measured experimentally. We note that since our external forcing P (or H) is in the wetting fluid, it is of opposite sign to the capillary pressure; thus, in our representation drainage occurs at lower P than imbibition, in contrast with the conventional capillary pressure curve.
Another interesting feature in Figure 4 is the behavior of the inner cycle. While the numerical simulations display perfect RPM (as proved in Holtzman et al. (2020)), the experimental inner cycle does not rejoin the primary drainage curve at the same point exactly. This is due to the larger steps of driving pressure (P or H) in the experiments ( Figure S1 in Supporting Information S1), deviating from the infinitesimal perturbations required in the quasistatic case, that could be realized numerically but not experimentally; this was confirmed by simulating coarser H increments.

Discussion
Our results indicate that the interfacial energy dissipated in Haines jumps accounts for most of the energy dissipated in slowly driven systems (e.g., as acoustic emissions ), in agreement with observations from experiments driven at constant rate (Berg et al., 2013;Hu et al., 2018;Måløy et al., 2021). For the experiment presented in this paper, our quasistatic approach accounts for ∼80% of the energy dissipated during fluid displacements ( Figure 4). We argue that the remaining ∼20% discrepancy originates from viscous dissipation caused by a finite velocity of the wetting fluid in experiments versus the zero velocity considered in the model. A nonzero velocity is an inherent feature of Haines jumps (rheons), regardless of how slow the interface is driven (Berg et al., 2013). Moreover, the corresponding viscous pressure drop depends on the sign of the front velocity, and thus plays opposite roles in drainage and imbibition. A small asymmetry between imbibition and drainage was already observed in experiments of slow displacements through a localized constriction (single defect) ; this asymmetry accumulates and intensifies in a disordered medium composed of multiple defects where deformations take place cooperatively as Haines jumps. The role of viscous dissipation in imbibition-drainage cycles at finite flow rate, and its impact on the property of RPM is the subject of ongoing research.
Another source of disparity between our model and quasistatically-driven experiments is the long relaxation timescale required for the interface to attain mechanical equilibrium at a new Jurin's height after each small pressure step (Clotet et al., 2012;Schlüter et al., 2017). In fact, the interface may continue to experience fluctuations even after reaching the Jurin's height (Lago & Araujo, 2001;Shikhmurzaev & Sprittles, 2012). Hence, the seemingly reversible, non-dissipative displacements (isons) also contribute to viscous dissipation, though to a lesser extent than Haines jumps (rheons). Additional mechanisms that can cause hysteresis and dissipation such as dynamic wetting, snap-off, and fluid trapping (Bonn et al., 2009;Giacomello et al., 2016; are not considered in our model, which is characterized by a single continuous, univalued interface. In conclusion, we have derived a thermodynamic framework that allows the quantification of energy dissipated via Haines jumps in capillary pressure-saturation trajectories at continuum scale from the microscopic mechanisms of surface tension and capillarity. The analysis presented here applies in fact to all modeling approaches  Helland et al., 2021). Return-point memory appears therefore as a very useful property to extend classical equilibrium thermodynamic principles to nonequilibrium systems driven through metastable equilibria. The present approach provides in this way a means of upscaling fluid-fluid displacements in disordered media, and a generic method of potential interest in quasistatically-driven disordered systems.

Data Availability Statement
All data used to generate the figures and conclusions in the paper can be downloaded from: https://dx.doi. org/10.6084/m9.figshare.22623262.