The origin of hysteresis and memory of two-phase flow in disordered media

Cyclic fluid-fluid displacements in disordered media feature hysteresis, multivaluedness, and memory properties in the pressure-saturation relationship. Quantitative understanding of the underlying pore-scale mechanisms and their extrapolation across scales constitutes a major challenge. Here we find that the capillary action of a single constriction in the fluid passage contains the key features of hysteresis. This insight forms the building block for an ab initio model that provides the quantitative link between the microscopic capillary physics, spatially-extended collective events (Haines jumps) and large-scale hysteresis. The mechanisms identified here apply to a broad range of problems in hydrology, geophysics and engineering. A wide variety of processes in nature and engineering ranging from rainfall and evaporation in soils to flow reversals in enhanced oil recovery and CO2 geosequestration exhibit hysteresis, multivaluedness, and memory in the imbibition/drainage cycles. Here, the authors investigate the origins of these phenomena by developing a model based on first principles that, by accounting for the impact of pore space microstructure, provides a link between microscopic capillary heterogeneity and large-scale hysteresis in pressure-saturation curves.

T wo-phase fluid displacements often occur in cycles, alternating between displacement of the less wetting fluid by the more wetting one (imbibition or wetting) and vice versa (drainage or dewetting). Examples include rain and evaporation in soils 1 , and flow reversal in enhanced oil recovery 2 or in CO 2 geosequestration 3 . Remarkably, drainage-imbibition cycles exhibit significant hysteresis. Pressure-saturation (PS) trajectories do not retrace the same path when the flow direction is reversed: the pressure required to achieve a given fluid saturation in drainage differs from that in imbibition 1 . Challenging our understanding of PS hysteresis at the continuum ("Darcy") scale is the multi-scale nature of two-fluid displacements in heterogeneous media: the behavior at the continuum scale is governed by micro-scale interface pinning and abrupt depinning events (instabilities or jumps) 4 . Due to interfacial tension, local instabilities become spatially correlated, giving rise to large, collective displacements involving large portions of the front (avalanches or Haines jumps), which appear to be a key mechanism in the origin of irreversibility, or hysteresis, in porous media flow [5][6][7] .
A remarkable feature of drainage-imbibition hysteresis is that it persists when displacements are carried out exceedingly slowly (quasistatically). Quasistatic hysteresis is common to a variety of disparate systems, from ferromagnetic materials to shapememory alloys, which have a rugged free-energy landscape with multiple local minima separated by large energy barriers. The minima correspond to long-lived metastable configurations in which, in the absence of fluctuations, the system remains trapped unless an external field drives it sequentially from one metastable configuration to the next. Quasistatic hysteresis is the macroscale manifestation of micro-scale energy-dissipative events along this evolution. It is also associated with history dependence, or memory, which reflects the inaccessibility of alternative states (energy minima) due to large energy barriers 8 . The rugged energy landscape results from structural disorder-either quenched-in or dynamically generated. In porous and fractured media it arises from the heterogeneous nature of the solid matrix.
Classical approaches to quasistatic hysteresis is "compartment" models, in which the macroscale response is the cumulative output of multiple independent hysteretic units called hysterons 7,[9][10][11][12][13] . Preisach and Everett models in ferromagnetism and porous media flow, respectively, belong to this general family. Compartment models can be tailored to reproduce the hysteretic behavior of a particular system, by deriving the suitable statistical distribution of hysterons from a complete set of first-order reversal PS trajectories 8,13 . However, they remain phenomenological in the sense that hysterons are introduced ad hoc, and the model parameters are not related to system properties such as permeability, wettability, and viscosity. Furthermore, compartment models consider all units of a given threshold to be invaded once the macroscopic forcing is exceeded, with no consideration of the spatial location of units with respect to the interface. Recently, models which include topological information in the form of Minkowski functionals were shown to capture multivaluedness in the PS relationship; however, at an appreciable computational cost [14][15][16] .
Here, we propose a model of quasistatic hysteresis that represents a significant improvement over phenomenological approaches such as compartment models. Our model applies to quasistatic two-phase displacements in a Hele-Shaw cell with random gap spacing. It relies on an accurate pressure balance equation for the fluid-fluid interface which simultaneously accounts for the competing effects of interfacial tension, hydrostatic pressure, and micro-scale capillary instabilities. In this way, it naturally adopts a form that resembles a fluctuationless random-field Ising model (RFIM) 17,18 , though written in terms of continuous variables. In contrast with phenomenological approaches, all parameters in our model have a clear, identifiable physical meaning, thus making it accessible eventually to detailed experimental scrutiny. The model captures interactions among neighboring heterogeneities through the in-plane interfacial tension, as well as the resulting lateral correlations that give rise to collective capillary rearrangements of the fluid-fluid interface between successive metastable equilibrium configurations (Haines jumps). At large scale, these collective events produce jerky PS curves featuring drainage-imbibition hysteresis cycles and returnpoint memory of internal PS trajectories, for which there is ample experimental evidence 5,19,20 .

Results
Interface model: theory and numerical implementation. To derive an interface model for quasistatic imbibition and drainage, we consider the following disordered medium: A thin open fracture with spatially variable aperture. A physical realization of this medium is an imperfect Hele-Shaw cell where localized gapthickness variations are produced by randomly placed heterogeneities or "defects" (constrictions and expansions; Fig. 1). Though the model applies generally to any two immiscible fluids of different density and viscosity, it is formulated below for a viscous liquid and a gas as wetting and nonwetting fluids. This enables comparison with the experiments in ref. 21 for an isolated defect.
The variability in gap thickness induces variability in capillary pressure, p c , through the change in the out-of-plane curvature of the meniscus. By contrast the interfacial tension induced by the in-plane curvature produces a restoring force through a localized Young-Laplace term that approximates curvature by the Laplacian of the interface height h(x). As explained in detail in Supplementary Note 1, the balance between these two pressure contributions and hydrostatic and driving pressures leads to the following stochastic differential equation for the equilibrium configurations of the interface height h(x) at a given external pressure head H (imposed pressure at the cell's inlet): where γ is the fluid-fluid surface tension, g the gravitational acceleration, and ρ the wetting fluid density. Here, we neglect gravitational effects in the nonwetting fluid. For a static contact angle θ the local capillary pressure is p c ðx; yÞ ¼ 2γ cos θ=bðx; yÞ, where b(x, y) = b 0 − δb(x, y) is the local gap thickness (b 0 is the unperturbed gap thickness, and δb > 0 and δb < 0 correspond to constrictions and expansions, respectively). The cell is tilted (by an angle α) to prevent viscous fingering in drainage, which would interfere with the effect of disorder. For compactness of notation we define an effective gravity g e ¼ g sin α, and the external pressure P = ρgH (Fig. 1). The Hamiltonian corresponding to Eq.
(1) is given by dy ρg e y À P À p c ðx; yÞ Â Ã : We define a local effective field p e (x) as the negative functional derivative of H with respect to h(x), p e ðxÞ ¼ ÀδH=δhðxÞ.
Equation (1) for the interface height is obtained by setting p e (x) = 0, such that h(x) minimizes the Hamiltonian H, see Supplementary Note 2. We note that the pressure balance equation (1) is a specific form of the quenched Edwards-Wilkinson (qEW) equation 22,23 in which the driving force is given by −ρg e h(x) + P. Upon a small change in P, the configuration of the interface changes until this driving force compensates the two pressure terms of capillary origin. At each step (P value), the average equilibrium position of the interface corresponds to a new Jurin height-where capillary rise is restrained by gravity. Each new Jurin height is then an absorbing state of the dynamics 24 . The qEW equation plays a salient role in the context of elastic manifolds in random media, as the representative model for a whole universality class of critical non-equilibrium interfacial problems [25][26][27] . Eq. (1) can be identified also with a first-order approximation in dh/dx of the fluctuationless RFIM for fluid fronts in disordered media, in the quasistatic limit [28][29][30][31] . These two models, the qEW equation and the fluctuationless RFIM, have been subjected to intense scrutiny regarding interfacial kinetic roughening and avalanche dynamics, and typically formulated in terms of effective parameters. Here by contrast all parameters in Eq. (1) have direct physical meaning. Moreover, to the best of our knowledge, this is the first time that an interfacial equation like Eq. (1) is used to investigate hysteresis and memory of interfacial configurations.
The numerical implementation of our model uses an iterative algorithm inspired by the method used for simulating the fluctuationless RFIM. We discretize the horizontal direction equidistantly, x n = nΔx. For a given external pressure P, we use synchronous updating to resolve the equilibrium configuration, corresponding to the vanishing local effective field p e (x) = 0: the local height h(x n ) is increased in imbibition by δh for all sites n which are out-of-equilibrium, p e (x n ) > 0. This procedure is repeated until p e (x n ) is everywhere smaller than a tolerance ϵ. During drainage the criterion p e (x n ) < 0 decreases the local interface height by δh until p e (x n ) > −ϵ. Details of the numerical implementation are given in "Methods".
Hysteresis across a single mesa defect. We begin by considering the case of a single mesa-shaped defect which locally reduces the gap thickness. This case has been studied theoretically and experimentally in ref. 21 as exponent of a fundamental mechanism for hysteresis that does not depend on disorder. Since it admits also an analytical solution (Supplementary Note 3) it serves here as a benchmark for our numerical implementation of Eq. (1). Figure 2a, b shows the simulated interface positions at various external pressures during imbibition and drainage, highlighting three specific external forcing values, H (marked 1-3; see Supplementary Movie 1 for the complete sequence). The interface is flat both before touching and after leaving the defect. During imbibition, the interface experiences an instability when it first touches the defect: the sudden increase of the out-of-plane curvature pulls the central part of the interface into the defect and forces it to deform abruptly. The abrupt change of the front upon touching the defect corresponds to the passage from a metastable configuration to a stable one when the former reaches its limit of metastability 21 . This instability manifests in an abrupt increase of saturation (relative volume fraction of the wetting phase, S w ) in the PS relationship (Fig. 2c). The deformed interface then maintains the same shape until it reaches the edge of the defect, where it remains pinned until it flattens and leaves the defect. This motion is fully reversible because as the interface exits the defect only one stable configuration exists 21 . During drainage, the meniscus is pinned upon contact with the defect until it recovers the fully deformed shape, and then advances keeping this shape (Fig. 2b). Upon reaching the edge of the defect it snaps irreversibly back to a flat interface. The jump in drainage occurs at a lower external pressure than in imbibition, giving rise to the elementary hysteresis cycle in the PS relation shown in Fig. 2c. The simulated interfaces for the three highlighted H values are compared in Fig. 2d, e (imbibition and drainage, respectively) with experiments from ref. 21 for identical system geometry and fluid properties, and with the analytical solution given in the "Methods". The numerical and analytical solutions are in perfect agreement, predicting the experimental data with no fitting parameters.
Hysteresis in disordered media. We now turn to disordered media, examining the collective effect of the interactions between multiple microscopic events associated with spatial heterogeneity (multiple defects) and its impact on the macroscopic PS relationship. Interfacial tension introduces spatial correlations through the in-plane interfacial curvature, such that capillaryinduced pinning or depinning of the front in one location affect their occurrence in other locations. This propagates in a cascade process that eventually modifies the configuration of the entire interface. Competition between interfacial tension and effective gravity (Eq. (1)) leads to a lateral correlation length ' g ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi γ=ðρg e Þ p which extends beyond the typical size of gap-thickness fluctuations 21 . This explains how strongly nonlinear behavior emerges from an equation that is derived for |dh/dx| ≪ 1.
The disorder introduced by randomly placed defects of variable thickness δb gives rise to a distribution of p c values. We separate p c into p 0 c þ δp c , where p 0 c ¼ 2γ cos θ=b 0 and δp c ¼ p 0 c ðδb=b 0 Þ=ð1 À δb=b 0 Þ. We model δp c as a random variable that takes on the value 0 with probability 1−ψ and a value drawn from a Gaussian distribution of zero mean and standard deviation σ 0 with probability ψ (where ψ represents the relative fraction of defects in the cell). This creates a rugged energy landscape which is manifested by spatially rough interfaces and Haines jumps. The latter are abrupt, notable changes to the interface configuration and therefore to the wetting-phase saturation, resulting in jerky PS behavior ( Fig. 3; see Supplementary Movie 2 for the dynamics).
The average height of the consecutive equilibrium configurations reached by the interface is given by the corresponding Jurin height. Hence, the average slope of the PS retention curves is dictated by the imposed effective gravity, whose relative importance in comparison with capillary forces is measured by a Bond number, Bo. In the limit of very small Bo (nearly horizontal cell) the dominance of capillarity allows for larger avalanches, and a larger change in saturation in response to a given change in external forcing, so that the PS curves become nearly horizontal, as observed experimentally 32 . For further analysis of the impact of Bo see Supplementary Note 5. In standard retention experiments, where the setup limits the sample size so that gravity becomes negligible relative to capillary forces (small Bo) the PS curves are nearly horizontal 1 .
A remarkable property of the model is its ability to capture the memory properties of PS trajectories. An internal hysteresis cycle rejoins the primary PS trajectory exactly at the point where the internal cycle was initiated ("return-point memory", RPM; Fig. 3a), with exactly the same interface shape. This is also true for cycles within internal cycles (not shown here), revealing that not only PS trajectories are infinitely multivalued, but also that their actual value depends on the previous sequence of return points in the current trajectory, which is stored in the interfacial configuration evolution (Fig. 3b). This RPM property is shared by a number of seemingly disparate slowly driven disordered systems with a rugged, complex free-energy landscape, from cellular automata 33 to ferromagnetic and shape-memory polycrystalline materials 34 . RPM has been observed experimentally as a robust feature in two-phase flows in porous media 5,19,20 . Our own preliminary experiments in a Hele-Shaw cell featuring randomly placed mesa-shaped constrictions confirm, within experimental accuracy, the numerically observed behaviors, also in terms of the RPM property. In contrast to the fluctuationless RFIM, our model involves continuous variables h(x) interacting through a Laplacian term. Nevertheless, the RPM property can be rigorously proved by analogous arguments (see Supplementary Note 4), bringing new insights about the interfacial dynamics of two-phase flows in disordered media implicit in the energy minimization leading to Eq. (1).
Hysteresis depends on the underlying microstructural heterogeneity in a non-trivial manner, as a result of the interactions among defects associated with surface tension. We consider the impact of the width of the distribution of p c , σ p c , on the strength of the hysteresis cycles, ΔP. The latter is measured by the average pressure difference between drainage and imbibition. We find that increasing σ p c results in stronger hysteresis (increasing ΔP); an example is shown in Fig. 4a (each curve is the ensembleaveraged PS values over 15 disorder realizations; see Supplementary Note 6). We find that ΔP grows as ΔP $ ðσ p c Þ n with n ≈ 1.4. Interestingly, this scaling holds for systems with different types of defect distributions (Fig. 4b). It also holds for: (1) dichotomic distributions of δp c with various relative fraction of defects ψ; and (2) with defects placed on a regular grid, with disorder introduced only through variations in δp c (see Supplementary Note 7). The common features of these different scenarios are the existence of a finite length scale characterizing the geometry relative to the correlation length ℓ g , and a finite defect thickness characterizing the disorder (e.g., its variance). In the single defect case, the saturation jump and therefore the hysteresis strength is linearly proportional to the defect strength, δp c 21 . The nonlinear dependence of ΔP on σ p c in disordered media is a result of interactions among the heterogeneities due to lateral correlation caused by interfacial tension. Note that the length of influence ℓ g (here~10 mm) is larger than the typical distance between the defects (~0.1 mm for ψ = 0.35 occupancy). ΔP is a measure of the energy barrier that needs to be overcome when alternating from imbibition to drainage. Our results point to a nonlinear dependence of this energy barrier on the width of the disorder distribution.

Discussion
Hysteresis is a consequence of energy dissipation; in our quasistatic approach the dissipation occurs through the irreversible release of elastic energy in Haines jumps 35 . A preliminary estimation based on Eq. (2) demonstrates that a single avalanche can dissipate a substantial amount of the energy stored in the deformed interface, in line with recent experimental   36 . Furthermore, since the RPM property ensures that the two-fluid interface recovers its original configuration in a closed imbibition-drainage cycle, it follows that the total energy dissipated in the cycle equals the total work made by the driving pressure in the cycle, i.e. the area enclosed in the PS cycle 37 . Thus, our results provide direct evidence that the dissipated energy increases with microstructural heterogeneity (Fig. 4).
Over the last century, fluid-fluid interface pinning via localized forces has been studied extensively, suggesting their intimate link to the measured (macroscopic) hysteresis. Yet, how this qualitative knowledge could be cast into a quantitative, physics-based predictive model remained elusive. Our work achieves that by deriving ab initio a model which relies on physical and measurable parameters only, allowing its rigorous validation against experimental and analytical data. Accounting for the impact of the microstructure, our model presents a significant stride from compartment models which consider disorder statistically, such that for a given external pressure all pores which are thermodynamically available in terms of capillary thresholds are filled regardless of their location relative to the interface 7,38 . In contrast, our model includes the essential distinction between availability and accessibility. We anticipate that the fundamental disorder and capillary mechanisms uncovered in this article can be transferred to more complex porous and fractured media. This new understanding will provide a much needed step change in our ability to design and monitor energy production, water resources and environmental engineering operations.

Methods
Numerical model implementation. The numerical solution of Eq. (1) for a given configuration of p c (x, y) is based on the minimization of the Hamiltonian, Eq. (2). Discretizing in space in the x direction, x k = kΔx, provides the following effective field (pressure imbalance) p e (x): where we set h(x k ) = h k . In imbibition (increasing h), after incrementing the external forcing P, the effective field is evaluated at all positions x k . At all the positions which are out-of-equilibrium, p e (x k ) > ϵ, the interface height h(x k ) is simultaneously updated-incremented by a constant δh. Note that in imbibition (drainage) p e is locally increased (decreased) when touching a positive defect, δb > 0, and decreased (increased) when touching a negative defect, δb < 0.
Advancing the interface at the location of the defect decreases p e (x k ). Due to surface tension (first term in Eq. (3)), this update affects neighboring interface positions, perturbing them away from equilibrium. Thus, when p e is again evaluated at all x k for the new interface configuration, the heights h(x k ) are updated if p e (x k ) > 0. These iterative search is repeated until equilibrium is met at all x k , p e (x k ) < ϵ. For drainage, the procedure is identical however with opposite signs: The interface height is decreased synchronously by δh at every out-of-equilibrium position x k , repeatedly until equilibrated everywhere, p e (x k ) > −ϵ. A pseudocode of the algorithm is provided as Supplementary Methods. The discretization Δx, the accuracy ϵ, and the height increment δh need to satisfy the following criteria. First, we note that very "weak" defects (small |δb|) for which |p c (x k , y)| < ϵ are invisible to this algorithm. Furthermore, one needs to ensure that the iterative increments in the gravitational potential energy and surface tension terms are sufficiently small, ρg e |δh| ≪ ϵ and γ|δh|/Δx 2 < ϵ, respectively. This means that |δh| ≪ ϵ/ρg e and Δx > (γ|δh|/ϵ) 1/2 . At the same time, Δx must be smaller than in order to resolve the interface curvature. Typically though, Δx will be smaller than the minimum defect width. If these conditions are not fulfilled, the effective field may accumulate negative (positive) values in imbibition (drainage). This in turn can lead to the following artifact: Upon reverting the flow direction and thus the equilibrium condition-from p e (x) < ϵ (imbibition) to p e (x) > −ϵ (drainage) and vice versa-the interface becomes immediately unstable and advances without changing the external forcing P. In all the simulations presented in this paper, we used steps of ΔH = 0.05 mm, numerical resolution of Δx = 0.1 mm and δh = 5 × 10−6 mm, and accuracy of ϵ = 0.01 Pa.
Experimental and numerical data for a single mesa defect. We perform an imbibition-drainage cycle by forcing oil in and out of a dry (air-filled) Hele-Shaw cell with a single mesa defect: a rectangular elevation of the bottom plate, locally  reducing the gap thickness b 0 by δb (Fig. 1). The cell is tilted at an angle α to avoid fingering instabilities. Interfacial instability in the form of viscous fingering is controlled by the sign of Bo−Ca, where Ca is the ratio of viscous to capillary forces 39 . Tilting the cell enforces Bo−Ca > 0, thus preventing fingering. The oil is forced in and out of the cell by raising and lowering an oil reservoir connected to the cell inlet (y = 0), which sets the pressure head H. The experiment in Fig. 2 is simulated with our model using the same conditions: w = 3 mm, α ¼ Analytical solution for a single mesa defect. In the following, we present analytical solutions for the interface height for a single mesa defect. For details of the derivations see Supplementary Note 3.
We consider a rectangular positive defect centered in x = 0, of width w and length ℓ (extending from y = d to y = d + ℓ). The defect strength is denoted by a 0 ¼ p 0 c δb=ðb 0 À δbÞ. For convenience, we define H e ¼ ðH þ p 0 c =ρgÞ= sin α and k ¼ ρg e =½1 À expðÀw=2' g Þ, where ' g ¼ ðγ=ρg e Þ 1=2 . During imbibition, while H e < d the interface is flat and h(x) = H e . After the meniscus touches the defect (H e ≥ d), capillarity sucks the interface into the defect and deforms it such that Equation (7) indicates that w e would increase with decreasing H e , which is not possible because it cannot exceed the physical width w. Therefore, the interface will snap-off at H c .

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The simulations used to produce the findings of this study are described in the pseudocode (Supplementary Methods), and the code is available from the corresponding author upon reasonable request. and defect width w ∈ [3,9,27] mm (the different symbols highlight the latter). The data from the 36 cases is collapsed by its presentation in a nondimensional form (some cases overlap), suggested by Eq. (6). The slope of the straight line is equal to cos θ (here cos θ ¼ 1). The excellent agreement between the simulated results (symbols) and the analytical, theoretical prediction (line) demonstrates our model accurate numerical implementation and convergence. It also predicts the corresponding experimental data, cf. Fig. 3 in ref. 21 .