#  Emergence of dissipation and hysteresis from interactions among reversible,
non-dissipative units: The case of fluid-fluid interfaces

Ran Holtzman<EMAIL_ADDRESS>Centre for Fluid and Complex
Systems, Coventry University, Coventry, United Kingdom  Marco Dentz  Institute
of Environmental Assessment and Water Research (IDAEA), Spanish National
Research Council (CSIC), Barcelona, Spain  Marcel Moura  PoreLab, The Njord
Centre, Department of Physics, University of Oslo, Oslo, Norway  Mykyta
Chubynsky  Centre for Fluid and Complex Systems, Coventry University,
Coventry, United Kingdom  Ramon Planet  Departament de Física de la Matèria
Condensada, Universitat de Barcelona, Barcelona, Spain, and University of
Barcelona Institute of Complex Systems, Barcelona, Spain  Jordi Ortín
Departament de Física de la Matèria Condensada, Universitat de Barcelona,
Barcelona, Spain, and University of Barcelona Institute of Complex Systems,
Barcelona, Spain

###### Abstract

We examine the nonequilibrium nature of two-phase fluid displacements in a
quasi-two-dimensional medium (a model open fracture), in the presence of
localized constrictions (“defects”), from a theoretical and numerical
standpoint. Our analysis predicts the capillary energy dissipated in abrupt
interfacial displacements (jumps) across defects, and relates it to the
corresponding hysteresis cycle, e.g. in pressure-saturation. We distinguish
between “weak” (reversible interface displacement, exhibiting no hysteresis
and dissipation) and “strong” (irreversible) defects. We expose the emergence
of dissipation and irreversibility caused by spatial interactions, mediated by
interfacial tension, among otherwise weak defects. We exemplify this
cooperative behavior for a pair of weak defects and establish a critical
separation distance, analytically and numerically, verified by a proof-of-
concept experiment.

##  I Introduction

Path-dependency (hysteresis) in pressure-saturation relationships during
imbibition and drainage in two-phase displacements in porous media occurs in
wide variety of natural and engineered processes, for instance soil moisture
and geoenergy [1, 2]. This hysteresis is largely due to individual and
cooperative capillary instabilities, known as Haines jumps [3], which are
inherently related to energy dissipation [4, 5]. The ubiquity of these
phenomena in various applications and the intriguing underlying physics
motivated extensive experimental, numerical and theoretical studies [5, 6, 7,
8, 9, 10, 11, 12]. The intrinsic complexity of porous media and non-local
pore-scale interactions makes the quantitative understanding of the precise
mechanisms that lead to these phenomena a challenging task. Modeling
approaches that account for the metastability of two-fluid configurations and
for hysteresis are often based on the aggregation of individual hysteretic
units (hysterons), in the so-called compartment models (e.g. [6, 7]).
Recently, based on quantitative insights from systematic studies in simplified
model systems that allow to isolate individual features causing hysteresis
[13, 14], a novel approach that does not rely on the concept of hysterons was
suggested in [15].

However, the possibility that non-hysteretic units interact cooperatively to
give rise to hysteresis and dissipation remains unexplored, despite evidences
from paradigmatic models of collective phenomena such as the Random Field
Ising Model (RFIM) [16, 17]. Here we address this open question by studying
the passage of a two-phase interface through localized perturbations in the
capillary pressure, representing single topographic defects (gap thickness
constrictions) in a 2-D medium (Hele-Shaw cell).

Two-phase fluid flow in Hele-Shaw cells—quasi-two-dimensional setups formed by
two closely spaced parallel plates [18, 19, 20]—shares some salient features
with two-phase flows in porous and fractured media. In both, the bulk behavior
of the flow of fluids in the viscous (Stokes) regime can be described by
Darcy’s law [21], and the stability of fluid-fluid interfaces is controlled by
viscosity and density contrasts between the fluids [22]. A richer, more
realistic model system is the “imperfect” Hele-Shaw cell, featuring gap-
thickness constrictions and expansions, in which the interfaces are subjected
to capillary instabilities akin to those observed in disordered media [13,
18]. Imperfect Hele-Shaw cells therefore stand out as an ideal playground to
study two-phase displacements in disordered media.

Using this model system, we study the notion of weak (reversible i.e. non-
dissipative and non-hysteretic) and strong (dissipative and hysteretic)
defects. We first classify single topographic defects as weak or strong, and
quantify the associated energy dissipation and hysteresis. We then
demonstrate, analytically, numerically, and experimentally, the emergence of
dissipation and hysteresis due to spatial interactions (surface tension) among
individually weak defects. We show that a pair of weak defects that is non-
dissipative if separated far apart becomes strong (dissipative) when brought
close enough together.

The paper is organized as follows. Section II provides background on
quasistatic pressure-driven displacements in imperfect Hele-Shaw cells,
including the mathematical statements of pressure equilibrium and energy
dissipation. Section III presents the solution for the interfacial
configurations across single defects (capillary pressure distortions). We show
the functional form of the interfacial shapes, which also reveal subtle
morphological differences between imbibition and drainage, and distinct
between weak (reversible, non-dissipative) and strong (irreversible,
dissipative) defects. Section IV uses the above to compute an energy balance
from which we establish the energy dissipated during jumps. In Section V we
consider cooperative effects in the presence of multiple defects, with a
nonintuitive result—energy dissipation due to spatially-correlated
interactions through interfacial tension for a pair of weak defects brought
sufficiently close together. A theoretical analysis provides a critical value
for the lateral separation between two weak defects that makes them
collectively dissipative. This prediction, derived in the limit of very narrow
defects, is validated against numerical simulations, and qualitatively
verified experimentally. Finally, in Sec. VI the original results of the work
are briefly recalled and the main conclusions are drawn. Specific details of
the calculations, numerical simulations and experiments are provided in the
Appendices.

##  II Background

###  II.1 Model system: Imperfect Hele-Shaw cell

Quasistatic pressure-driven displacements in imperfect Hele-Shaw cells have
shown strongly nonlinear behavior and Haines jumps at pore and multipore
scale, together with lack of reversibility between opposite displacement
directions (imbibition and drainage) at continuum scale, hysteresis in the
applied pressure vs wetting-phase saturation (PS) trajectories, and the
return-point memory (RPM) property of closed partial cycles that is ubiquitous
in porous media flows [14, 15, 23]. Below, we describe the theory allowing to
compute the interfacial configurations and energy dissipation for the general
case of an imperfect Hele-Shaw cell with multiple defects (representing
disordered media). In this paper we apply these concepts to displacements
through (i) isolated defects (local constrictions), and (ii) a pair of
defects, exposing their interaction.

###  II.2 Establishing the equilibrium configurations

The pressure balance of a two-phase interface invading an imperfect Hele-Shaw
cell can be derived in quasistatic conditions [15, 14]. Without loss of
generality, we assume in the following that the less wetting fluid is low-
viscosity (e.g. air) and the more wetting fluid is viscous (e.g. liquid such
as silicone oil). We consider that (i) the gap thickness in the cell changes
in space (giving rise to an extended domain of connected constrictions and
expansions); (ii) the fluids are immiscible, and displacements are driven by
the change of the imposed pressure $P$ at one end of the cell; and (iii) the
cell is tilted in the direction that prevents the formation of viscous fingers
[22]. With these conditions, the linearized pressure balance takes the form

$\gamma\frac{d^{2}h(x)}{dx^{2}}-\rho g_{e}h(x)+P+p_{c}\left[x,h(x)\right]=0,$
(1)

where $h(x)$ is the equilibrium interface position at $x$, $\gamma$ the oil-
air surface tension, $\rho$ the oil density, $g_{e}$ is the effective gravity
(which in a physical implementation could be changed by tilting the cell,
$g_{e}=g\sin{\alpha}$, where $g$ is the gravitational acceleration and
$\alpha$ the inclination angle from the horizontal). Here $p_{c}(x,y)$ is the
perturbation in out-of-plane capillary pressure, determined by the variations
in thickness ($z$). Here, the direction of the fluid advancement is parallel
to the $y$-axis. For simplicity, we do not account for the minute variations
in hydrostatic pressure of non-wetting fluid relative to the liquid pressure
of the wetting fluid. Additionally, in our quasi-static model, the pressure
changes associated with the viscosity of the two fluids are also neglected.
The first and last terms in Eq. (1) account for the linearized in-plane
component of the Young-Laplace pressure jump across the interface at each site
$x$ (for comparison with the exact nonlinear term see [23]), and the out-of-
plane component arising from the presence of expansions and constrictions in
the cell, respectively. The role of the two terms is different: the out-of-
plane component is responsible for the interface deformation whereas the in-
plane component is a restoring force resisting the deformation. In the
quasistatic limit (zero driving rate) displacements are driven by minute
changes of $P$, separated by long time intervals required for reaching a new
mechanical equilibrium $h(x)$.

The equilibrium configurations $h(x)$ could also be derived from minimizing
the Hamiltonian

$\mathcal{H}=\int_{-\infty}^{\infty}dx\left(\frac{\gamma}{2}\left[\frac{dh(x)}{dx}\right]^{2}+\int\limits_{0}^{h(x)}dy\left[\rho
g_{e}y-P-p_{c}(x,y)\right]\right),$ (2)

and the pressure imbalance $p_{e}(x)$ experienced by the interface at each
site $x$ is given by $p_{e}(x)=-{\delta\mathcal{H}}/{\delta h(x)}$. The
condition of mechanical equilibrium in Eq. (1) corresponds therefore to
setting $p_{e}(x)=0$ [15]. In the framework of this model, PS trajectories are
built from the sequence of equilibria. The passage from one equilibrium
configuration to the next can be of two kinds. In the first, the system
remains trapped in a _local_ energy minimum, where the small change of
external forcing $P$ causes a correspondingly small evolution of the wetting-
phase saturation $S_{w}$, resulting in a smooth PS trajectory. In the second,
an abrupt change of state ($S_{w}$) takes place at the new value of $P$, in a
Haines jump [3, 5]; this occurs when the change of $P$ suppresses the current
local energy minimum, and the system is forced to abruptly jump to a new
metastable equilibrium. Haines jumps are effectively instantaneous in the time
scale of change of the driving pressure, so that interfacial configurations
experience irreversible changes at punctuated values of $P$.

Numerically, this dynamics can be simulated iteratively by _synchronous_
updates of $h(x)$ in all unstable sites by a small amount in the direction
that reduces $|p_{e}(x)|$, stopping when all sites retain equilibrium (for
details of the numerical procedure, see [15]). This deterministic rule is akin
to the zero-temperature limit of the Glauber dynamics for RFIM [16, 24], which
considers energy barriers between consecutive equilibria that are much larger
than thermal fluctuations. The presence of the quenched disorder term
$p_{c}[x,h(x)]$ in Eqs. (1) and (2) defines a rugged free energy landscape, so
that for every applied pressure there are many different interfacial
configurations $h(x)$ that are local minimizers of $\cal{H}$. The synchronous
dynamics described above takes the current configuration to the closest
available metastable minimum in a deterministic manner dictated by the
quenched disorder. Ref. [15] proved that no parts of the interface recede
under this dynamics, and a no-passing rule [25] is obeyed such that a
configuration of larger (or equal) saturation compared to another will remain
so under a monotonous evolution of the driving pressure. As a result, the
original two-phase configuration is exactly recovered in any cyclic excursion
of the driving pressure, a property known as RPM (return-point memory),
ubiquitous to many athermal driven disordered systems [26, 27, 28, 29].

###  II.3 Energy dissipated between equilibrium configurations

The amount of energy dissipated can be obtained from the change in interfacial
energy due to fluid displacement and the mechanical work done by the applied
pressure $P$. For a small interface displacement $\delta h(x)$ the dissipated
energy is

$d\bar{}\hskip 1.00006pt\Psi=dU-d\bar{}\hskip 1.00006ptW,$ (3)

where $dU$ is the change in the internal energy, and

$d\bar{}\hskip 1.00006ptW=PdS=P\int_{-\infty}^{\infty}\delta h(x)dx$ (4)

is the work. Here the notation $d\Box$ is used for infinitesimal changes of
variables that are state functions (e.g. $S$), while $d\bar{}\hskip
1.00006pt\Box$ is for changes of variables that are not (e.g. $W$). In our 2-D
model, all energy units [Eq. (4) and throughout] are of energy per unit
length, i.e. the out-of-plane thickness. $U$ accounts for the capillary energy
of the front deformation and for the gravitational potential energy of the oil
phase (wetting fluid). We follow the convention that $d\bar{}\hskip
1.00006pt\Psi\leq 0$.

Noting that $\mathcal{H}=U-PS$, we get $\delta\mathcal{H}=d\bar{}\hskip
1.00006pt\Psi-S\,dP$. The fact that
${\mathcal{H}}={\mathcal{H}}\left[h(x),P\right]$ and $p_{e}=-\delta{\cal
H}/\delta h(x)$ allows writing

$\delta{\mathcal{H}}=-\int_{-\infty}^{\infty}dx\,p_{e}(x)\,\delta
h(x)-\int_{-\infty}^{\infty}dx\int_{0}^{h(x)}dy\,dP.$ (5)

From this, as well as from Eq. (2) and the expression ${\cal H}=U-PS$, we can
compute the internal energy as

$U=\int\limits_{-\infty}^{\infty}dx\left[\frac{\gamma}{2}\left(\frac{\partial
h}{\partial x}\right)^{2}+\int\limits_{0}^{h}dy\left(\rho
g_{e}y-p_{c}\right)\right]$ (6)

Finally, the energy lost between two equilibrium states, $t-1$ and $t$, is
found by integrating Eq. (3), which gives [12]

$\Psi^{t-1\to
t}=\left[U^{t}-U^{t-1}\right]-P^{t}\left[S_{w}^{t}-S_{w}^{t-1}\right].$ (7)

###  II.4 Energy dissipated per incremental change in interface configurations

Here, we present an alternative method for computing energy dissipation
between consecutive interfacial configurations. For continuous, reversible
displacements (isons), where the system stays in one local minimum of
$\mathcal{H}$, and a small change in $P$ leads to a small change in $h(x)$,
$p_{e}(x)=0$ by definition. This reduces Eq. (5) to
$\delta{\mathcal{H}}=-\int_{-\infty}^{\infty}dx\int_{0}^{h(x)}dy\,dP=-S\,dP,$
which proves that $d\bar{}\hskip 1.00006pt\Psi=0$, i.e., no dissipation. Since
there are no energy losses in this case, $d\bar{}\hskip 1.00006ptW=dU$, and
for a finite continuous displacement $W=\Delta U$, with $\Delta U$ computed
from Eq. (6).

In irreversible displacements (rheons) $P$ is constant ($dP=0$) and
$\delta\mathcal{H}=d\bar{}\hskip 1.00006pt\Psi$. This in turn provides the
energy dissipated for each elementary step $\delta h$ within a Haines jump,
using Eq. (5):

$d\bar{}\hskip
1.00006pt\Psi=\delta{\mathcal{H}}=-\int_{-\infty}^{\infty}dx\,p_{e}(x)\,\delta
h(x).$ (8)

To obtain the total energy dissipated in a given jump between equilibrium
states $t-1$ and $t$, $\Psi^{t-1\to t}$, one can integrate Eq. (8) over all
the intermediate nonequilibrium steps $k$,

$\Psi^{t-1\to t}=-\sum_{k}\int_{-\infty}^{\infty}dx\,p_{e}^{k}(x)\,\delta
h^{k}(x).$ (9)

Equations (8) and (9) present the energy dissipated as the sum of elementary
products of unbalanced pressure $p_{e}(x)$ times the corresponding
displacement $\delta h(x)$. This will be used below (Section III) to provide
an explanation for a nonintuitive observation: asymmetry between energy
dissipation in drainage and imbibition. Furthermore, Eqs. (8)–(9) show that
the energy dissipation within avalanches does not have to be proportional to
the corresponding change in saturation, as the values of $p_{e}(x)$ can be
different from site to site. This non-proportionality between changes in
saturation and dissipation was shown numerically in disordered media
constructed from defects of various strengths [12]. A related finding was
shown for quasistatically driven disordered ferromagnets [17]. An extreme
example of this non-proportionality arises in the limit in which dissipation
approaches zero even as the size of the interface jump remains finite, as we
will see below.

##  III Interface configurations: Single defect

In this Section, we formulate an analytical solution for the equilibrium
interface configuration for a single defect. When an incompressible wetting
fluid (e.g. oil) imbibes or drains quasistatically in a smooth Hele-Shaw cell
with a narrow gap of fixed width, filled with inviscid, non-wetting fluid
(e.g. air), the interface is morphologically stable and $h(x)=h_{0}$.
Modulations of the out-of-plane capillary pressure,
$p_{c}(x,y)=p_{c}^{0}+\delta p_{c}(x,y)$, cause the interface to deform [30,
14, 15]. Expressing the modulation as $\delta p_{c}(x,y)=\delta
p_{c}^{*}F(x,y)$, where $\delta p_{c}^{*}$ is the maximum value and $F(x,y)$
is the “normalized modulation”, we can rewrite Eq. (1) as [14]

$\gamma\frac{d^{2}\eta(x)}{dx^{2}}-\rho g_{e}\eta(x)+\delta
p_{c}^{*}F[x,\eta(x)+h_{0}]=0\,,$ (10)

where $\eta(x)=h(x)-h_{0}$ (Fig 1). Note that $h_{0}=(P+p_{c}^{0})/\rho
g_{e}$. In the following we consider the case of $\delta p_{c}^{*}>0$ ($\delta
p_{c}^{*}<0$ is analogous). We consider a defect with $F(x,y)$ nonzero within
a rectangle $-w/2<x<w/2$, $y_{d}<y<y_{d}+\ell$, where $w$ is the width of the
defect and $\ell$ its length. We do not set $y_{d}=0$ to keep our formulation
general such that it can be used for a disordered system with multiple
defects.

Figure 1:  Our 2-D model system representing an imperfect Hele-Shaw cell with
local perturbations in thickness of width $w$ and length $\ell$ (defects; in
brown), which alter the out-of-plane capillary pressure $p_{c}$. Once the
fluid-fluid interface (blue line) passes through a defect, the interface
deforms by $\eta(x)$, measured relative to the flanks (which are at $y=h_{0}$)
such that the interface height $h(x)=h_{0}+\eta(x)$ is nonuniform. We consider
both “mesa” defects (fixed $p_{c}$ within the defect) and regular defects
which include a sloping part of length $\ell_{s}$ (linear increase in $p_{c}$)
and a plateau (fixed $p_{c}$). We also consider both a single, isolated
defect, and a pair of interacting defects (separated by $d$); for simplicity,
the interface, as well as its descriptors [$h_{0}$, $\eta(x)$], and the defect
position $y_{d}$ are not shown for the defect pair.

The effective pressure field given by the left-hand side of Eq. (10) can be
split into two parts: $p_{e}(x)=p_{d}(x)+\delta p_{c}^{*}F[x,\eta(x)+h_{0}]$,
where $p_{d}$ accounts for the restoring force of the line and it is linear in
$\eta$ [14]. We obtain the equilibrium states, $p_{e}=0$, by equating

$p_{d}(x)=-\delta p_{c}^{*}F[x,\eta(x)+h_{0}].$ (11)

Here, we find two different scenarios depending on the number of possible
equilibrium configurations (roots) $\eta_{c}$ that fulfill Eq. (11). If for
all $h_{0}$ there is only one root, the defect does not lead to hysteresis or
dissipation and is termed “weak”. In contrast, for an hysteretic and
dissipative (“strong”) defect, for some interval of $h_{0}$ the solution gives
three equilibrium points: two stable equilibrium configurations, with the
largest and smallest $\eta$, while the intermediate one is unstable. We note
that for complex $p_{c}$ profiles it is possible to have more than three
roots; here we consider only linear variations in $p_{c}$.

Figure 2: Geometrical construction of the equilibrium positions of the
interfaces in the presence of a localized defect: (a) weak; (b) strong; and
(c) mesa. Blue lines correspond to the pressure field exerted by the defect.
Orange lines correspond to the restoring force (pressure) of the interface,
opposing the deformation. At equilibrium, $p_{d}=p_{c}$, and these states are
presented as black dots for stable equilibrium. Red dots correspond to
unstable states.

Figure 2 shows visual representations of these different scenarios. These 1-D
representations are exact for infinitely wide defects (when there is no $x$
dependence), and remain qualitatively valid for defects of arbitrary width. We
refer to panels (a–b) as “regular” defects, where the change of $p_{c}$ is
continuous along the $y$ axis. Panel (c) shows the special case of the “mesa”
defect [13], where the modulation $F[x,\eta(x)]$ is a rectangular function in
$y$, so that the change in $p_{c}$ is discontinuous (non-regular). The 2-D
interface shape when passing through a gap modulation can be derived by
realizing that Eq. (10) can be written as $\delta
p_{c}^{*}F[x,\eta(x)+h_{0}]=-\mathcal{L}\eta(x)$ with a linear differential
operator $\mathcal{L}$. The interface shape $h(x)$ is obtained using the
Green’s function formalism [14]

$\eta(x)=\int\limits_{-\infty}^{\infty}dx^{\prime}G_{0}(x-x^{\prime})\delta
p_{c}^{*}F[x^{\prime},h_{0}+\eta(x^{\prime})],$ (12)

where $G_{0}(x)$ satisfies

$\gamma\frac{d^{2}G_{0}(x)}{dx^{2}}-\rho g_{e}G_{0}(x)=\delta(x)\,,$ (13)

It takes the form [14]

$G_{0}(x)=\frac{\ell_{c}}{2\gamma}\exp(-|x|/\ell_{c}),$ (14)

with $\ell_{c}=\sqrt{\gamma/\left(\rho g_{e}\right)}$.

###  III.1 Interface configurations: Mesa defect

For the special case of a mesa defect,

$F(x,y)=\displaystyle{\rm Rect}\left(\frac{x}{w}\right){\rm
Rect}\left(\frac{y-(y_{d}+\ell/2)}{\ell}\right)\,,$ (15)

where ${\rm Rect}(x)$ is the rectangular function ${\rm
Rect}(x)=H(x+1/2)-H(x-1/2)$, $H(x)$ is the Heaviside function, $w$ is the
width and $\ell$ the length of the defect. When passing such defects, the
interface displacement is always hysteretic. The analytical solutions for the
interface shape when passing through a mesa defect, derived in [15], are
summarized below for completeness.

Assuming that for $|x|<w/2$, $h_{0}+\eta$ is between $y_{d}$ and $y_{d}+\ell$,
the integral in Eq. (12) can be evaluated analytically, which gives

$\eta(x)=\frac{\delta p_{c}^{*}}{\rho
g_{e}}\begin{cases}\exp(-|x|/\ell_{c})\sinh(w/2\ell_{c})\,,&|x|>w/2\\\
1-\exp(-w/2\ell_{c})\cosh(x/\ell_{c})\,,&|x|\leq w/2\end{cases}$ (16)

The maximum deformation occurs at the defect centre $x=0$, and equals

$\eta_{m}=\frac{\delta p_{c}^{*}}{\rho
g_{e}}\left[1-\exp(-w/2\ell_{c})\right]\,.$ (17)

During imbibition, the interface remains flat until it contacts the defect,
$h_{0}=y_{d}$, at which point it deforms abruptly to a shape given by Eq.
(16); here and elsewhere in this paper, we assume that the defect is
sufficiently long to accommodate the deformed interface, i.e.,
$\ell>\eta_{m}$. At the end of the defect, $y=y_{d}+\ell$, the interface gets
pinned and deformation decreases continuously to zero (flat) by the time
$h_{0}=y_{d}+\ell$, at which point the interface exits the defect. During
drainage, the initial part of the interface displacement is reversible
(identical to that in imbibition), i.e. it gets pinned and deforms until
$h_{0}=(y_{d}+\ell)-(\delta p_{c}^{*}/\rho
g_{e})\left[1-\exp(-w/2\ell_{c})\right]$. However, when $h_{0}$ reaches
$y_{d}$ (the point of interface jump in imbibition), the trivial solution
$\eta=0$ appears as a metastable solution, but the interface remains deformed
(thus, hysteresis is observed). The jump is delayed until $h_{0}$ is equal to
the external head $h_{c}$, given by

$h_{c}=y_{d}-\frac{\delta p_{c}^{*}}{2\rho
g_{e}}\left[1-\exp(-w/\ell_{c})\right]\,.$ (18)

At that point the interface passes through the corners of the defect. For
lower $h_{0}$, the effective width $w_{e}$ across which the defect is wet
would have to be smaller than $w$, and the new equilibrium position would be

$h_{0}=y_{d}-\frac{\delta p_{c}^{*}}{2\rho
g_{e}}\left[1-\exp(-w_{e}/\ell_{c})\right]>h_{c}.$ (19)

However, this is not possible [15] and the nontrivial solution disappears.

Note that for a very wide defect ($w\gg\ell_{c}$), immediately before the jump
during drainage the interface profile around $x=-w/2$ can be approximated as

$h(x)={y_{d}+}\frac{\delta p_{c}^{*}}{2\rho g_{e}}\begin{cases}\exp(\Delta
x/\ell_{c})-1\ ,&\Delta x<0,\\\ 1-\exp(-\Delta x/\ell_{c})\ ,&\Delta
x>0.\end{cases}$ (20)

where $\Delta x=x+w/2$ and $|\Delta x|\ll w$. Thus, the solution is
approximately symmetric with respect to the corner of the defect ($x=-w/2$,
$y=y_{d}$) as the center of symmetry [$h(-w/2+\Delta
x)-y_{d}\approx-\\{h(-w/2-\Delta x)-y_{d}\\}$]. Similarly for $x=w/2$ and the
corresponding corner. We will compare this to the case of a wide regular
defect with a slope below.

###  III.2 Interface configurations: Regular defect

We now consider a defect with a capillary pressure profile which varies in
$y$,

$F(x,y)=a_{1}(x)c(y),$ (21)

where

$a_{1}(x)={\rm Rect}\left(\frac{x}{w}\right)$ (22)

and

$c(y)=\frac{1}{\ell_{s}}(y-y_{d}){\rm
Rect}\left(\frac{y-y_{d}-\ell_{s}/2}{\ell_{s}}\right)+{\rm
Rect}\left(\frac{y-y_{d}-(\ell_{s}+\ell)/2}{\ell-\ell_{s}}\right).$ (23)

The defect has width $w$ and length $\ell$, with the capillary pressure
profile consisting of two parts: a slope (ramp) of length $\ell_{s}$ where
$p_{c}$ increases linearly in $y$ and a plateau of fixed $p_{c}$ (as in the
mesa defect, cf. Fig. 1). In the ramp, the pressure slope is

$\Pi=\frac{dp_{c}(0,y)}{dy}=\frac{\delta p_{c}^{*}}{\ell_{s}}.$ (24)

The interface deformation is calculated by substituting Eq. (21) in the
integral equation (12),

$\eta(x)=\delta
p_{c}^{*}\int\limits_{-\infty}^{\infty}dx^{\prime}G_{0}(x-x^{\prime})a_{1}(x^{\prime})c[h_{0}+\eta(x^{\prime})],$
(25)

where the propagator $G_{0}(x)$ is given by Eq. (14).

In imbibition, prior to touching the defect, $h_{0}<y_{d}$, the interface is
undeformed, $\eta(x)=0$. As the interface enters the defect, $h_{0}>y_{d}$, as
long as the interface deformation is small enough such that inside the defect
it remains within the ramp, Eq. (25) becomes

$\eta(x)=\Pi\int\limits_{-\infty}^{\infty}dx^{\prime}G_{0}(x-x^{\prime})a_{1}(x^{\prime})(h_{0}-y_{d})+\Pi\int\limits_{-\infty}^{\infty}dx^{\prime}G_{0}(x-x^{\prime})a_{1}(x^{\prime})\eta(x^{\prime}).$
(26)

We distinguish between two cases: the defect is _weak_ if a solution of this
equation exists, at least for small enough positive $h_{0}-y_{d}$; it then
changes continuously as $h_{0}$ increases. The defect is _strong_ if there are
no solutions for positive $h_{0}-y_{d}$; in that case, upon entering the
defect the interface experiences an abrupt jump that extends into the plateau.

To obtain a closed-form analytical solution, we consider below two limits:
infinitely narrow and infinitely wide defects. For these, we find the
threshold between weak and strong, and the interface shapes before and after
the jump in the strong regime. For the general case of a finite defect width,
the threshold calculation is provided in Appendix A, and the numerical methods
for establishing the interface shape are described in Appendix B.

The threshold between weak and strong defects in terms of the pressure slope,
$\Pi_{c}^{(1)}$, is independent of $\ell_{s}$ given $\ell_{c}$, $w$ and
$\gamma$. Dimensionality considerations provide

$\Pi_{c}^{(1)}=\frac{\gamma}{\ell_{c}^{2}}f\left(\frac{w}{\ell_{c}}\right),$
(27)

where $f$ is a dimensionless function. Given $\Pi/\Pi_{c}^{(1)}$ and
$w/\ell_{c}$, the interface deformations before the jump in imbibition,
$\eta_{\rm imb}$, and after the jump in drainage, $\eta_{\rm dra}$, are
proportional to $\ell_{s}$, thus,

$\eta_{\rm imb}(x)=\ell_{s}\phi_{\rm
imb}\left(\frac{x}{\ell_{c}};\frac{\Pi}{\Pi_{c}^{(1)}},\frac{w}{\ell_{c}}\right)$
(28)

and

$\eta_{\rm dra}(x)=\ell_{s}\phi_{\rm
dra}\left(\frac{x}{\ell_{c}};\frac{\Pi}{\Pi_{c}^{(1)}},\frac{w}{\ell_{c}}\right),$
(29)

where $\phi_{\rm imb}$ and $\phi_{\rm dra}$ are dimensionless functions.

####  III.2.1 Narrow defect limit ($w\ll\ell_{c}$)

Consider a narrow defect, $w\ll\ell_{c}$, located at $x=0$. We note that
$a_{1}(x)=w\delta_{w}(x)$ with $\delta_{w}(x)=w^{-1}{\rm Rect}(x/w)$, where in
the limit $w\to 0$, $\delta_{w}(x)$ approaches the Dirac delta $\delta(x)$.
Introducing this approximation in Eq. (26), we obtain

$\eta(x)=\Pi w(h_{0}-y_{d})G_{0}(x)+\Pi w\eta_{0}G_{0}(x),$ (30)

where $\eta_{0}=\eta(0)$. For $x=0$, this becomes

$\eta_{0}=\Pi w(h_{0}-y_{d})\beta_{0}+\Pi w\eta_{0}\beta_{0},$ (31)

with $\beta_{0}=G_{0}(0)=\ell_{c}/2\gamma$. This gives

$\eta_{0}=\frac{\Pi w\beta_{0}(h_{0}-y_{d})}{1-\Pi w\beta_{0}},$ (32)

so that we obtain the following solution:

$\eta(x)=G_{0}(x)\frac{\Pi w(h_{0}-y_{d})}{1-\Pi w\beta_{0}}.$ (33)

This is consistent with our assumption that the solution for the interface
configuration crosses the defect within the ramp when
$y_{d}<h_{0}+\eta_{0}<y_{d}+\ell_{s}$. This condition is satisfied for
sufficiently small and positive $h_{0}-y_{d}$, when the denominator of Eq.
(33) is positive, i.e.,

$\Pi<\Pi_{n}^{(1)}=\frac{1}{w\beta_{0}}=\frac{2\gamma}{w\ell_{c}}.$ (34)

In this case, the defect is weak. Here $\Pi_{n}^{(1)}$ is the narrow-defect
approximation for the critical slope in a single defect, where the general
threshold for a single defect, $\Pi_{c}^{(1)}$, is derived in Appendix A.
Conversely, for $\Pi>\Pi_{n}^{(1)}$ there is no solution crossing the defect
within the ramp for positive $h_{0}-y_{d}$, indicating that the interface
deforms abruptly beyond the sloping part, and the defect is strong. Note that
a solution crossing the defect within the ramp does exist for strong defects
when $h_{0}-y_{d}$ is small enough and negative; this is the unstable solution
marked by the red dot in Fig. 2(b).

For the case where the interface deforms beyond the ramp length and reaches
the plateau region, Eq. (25) gives

$\eta(x)=\delta p_{c}^{*}wG_{0}(x)=\delta
p_{c}^{*}\frac{\ell_{c}w}{2\gamma}\exp(-|x|/\ell_{c})=\frac{\Pi}{\Pi_{n}^{(1)}}\ell_{s}\exp(-|x|/\ell_{c}).$
(35)

In cases where the defect width cannot be neglected, we make the assumption
that the deformation is constant within the defect, providing

$\eta(x)=\frac{\Pi}{\Pi_{n}^{(1)}}\ell_{s}\begin{cases}\exp[-(|x|-w/2)/\ell_{c}],&|x|>w/2,\\\
1,&|x|<w/2.\end{cases}$ (36)

As expected, this coincides with the $w\ll\ell_{c}$ limit of the mesa defect
case, Eq. (16). The deformation is independent of $h_{0}$ as long as the
interface solution crosses the defect within the plateau. For a strong defect
($\Pi/\Pi_{n}^{(1)}>1$), this occurs already for $h_{0}=y_{d}$ [again,
considering a sufficiently long defect,
$\eta(0)=(\Pi/\Pi_{n}^{(1)})\ell_{s}<\ell$]. Thus, during imbibition, the
interface jumps from

$h=h_{\rm imb}^{-}=y_{d}$ (37)

to

$h=h_{\rm
imb}^{+}=y_{d}+\frac{\Pi}{\Pi_{n}^{(1)}}\ell_{s}\begin{cases}\exp[-(|x|-w/2)/\ell_{c}],&|x|>w/2,\\\
1,&|x|<w/2.\end{cases}$ (38)

In drainage, in the beginning the interface displacement is identical to the
mesa case, i.e., it pins at the edge of the defect, $y=y_{d}+\ell$, until the
deformation reaches the value given by Eq. (36), and then moves continuously.
The jump in drainage occurs when the interface reaches the boundary between
the plateau and the ramp (rather than the end of the defect, $y=y_{d}$, in the
mesa case), and the solution for the interface configuration, (36), meets the
unstable branch and ceases to exist. At this point,
$h_{0}={y_{d}+}\ell_{s}-\eta(0)={y_{d}}-\ell_{s}(\Pi/\Pi_{n}^{(1)}-1)$, and
the jump occurs between

$h_{\rm
dra}^{-}=\begin{cases}y_{d}-\ell_{s}\left\\{{\displaystyle\frac{\Pi}{\Pi_{n}^{(1)}}}[1-\exp(-(|x|-w/2)/\ell_{c})]-1\right\\},&|x|>w/2,\\\
y_{d}+\ell_{s},&|x|<w/2\end{cases}$ (39)

and

$h_{\rm dra}^{+}=y_{d}-\ell_{s}(\Pi/\Pi_{n}^{(1)}-1).$ (40)

These interfacial jumps at various defect strengths, computed numerically for
a sufficiently narrow defect $w/\ell_{c}=0.2$ (Appendix B), are illustrated in
Fig. 3a. As the transition between weak and strong is approached,
$\Pi/\Pi_{n}^{(1)}\to 1+$, the difference between $h_{\rm imb}^{+}$ and
$h_{\rm dra}^{-}$ vanishes (barely noticeable for $\Pi/\Pi_{n}^{(1)}=1.01$),
although the jump remains finite (noting it can vanish for other defect
profiles, not considered here).

Figure 3:  Interface profiles before and after jumps for a narrow
[$w/\ell_{c}=0.2$; (a)] and a wide [$w/\ell_{c}=20$; (b)] regular ramp
defects, for different capillary pressure slopes. The ramp part of the defect
is shown by a shading gradient, and the plateau part in uniform gray. The
boundary between these two regions is at $h/\ell_{c}=1$ and in panel (b) is
marked with a green dotted line. Dashed and solid lines represent the
undeformed (flat) and deformed configurations, i.e. before and after the jump
in imbibition, and the opposite in drainage.

####  III.2.2 Wide defects ($w\gg\ell_{c}$)

In case of a very wide defect, the width of the Green’s function $G_{0}$
[which is on the order of $\ell_{c}$; see Eq. (14)] is much smaller than the
width of the defect, and $G_{0}$ can be approximated by a delta function,

$G_{0}(x)\approx\frac{\ell_{c}^{2}}{\gamma}\delta(x),$ (41)

such that Eq. (12) becomes

$\eta(x)=\frac{\ell_{c}^{2}}{\gamma}\delta p_{c}^{*}F[x,h_{0}+\eta(x)].$ (42)

This approximation is valid far from lateral boundaries of the defect. Outside
the defect, $\eta(x)=0$. If the interface solution crosses the defect within
the ramp, then inside the defect we get from the above

$\eta=\frac{\ell_{c}^{2}\Pi}{\gamma}(h_{0}+\eta-y_{d}),$ (43)

and from this, or, equivalently from Eq. (26) with the same approximation for
$G_{0}$,

$\eta=\frac{(\ell_{c}^{2}\Pi/\gamma)(h_{0}-y_{d})}{1-\ell_{c}^{2}\Pi/\gamma}=\frac{\Pi(h_{0}-y_{d})}{\rho
g_{e}-\Pi}.$ (44)

This solution crosses the defect within the ramp for small enough positive
$h_{0}-y_{d}$ if $\Pi<\Pi_{w}^{(1)}=\rho g_{e}={\gamma}/{\ell_{c}^{2}}$. In
this case, the defect is weak. Here $\Pi_{w}^{(1)}$ is the wide-defect
approximation for the general solution for the threshold for a single defect,
$\Pi_{c}^{(1)}$, derived in Appendix A. Otherwise (for $\Pi>\Pi_{w}^{(1)}$),
during imbibition the interface jumps into the plateau upon touching it at
$y_{d}$. In that case, the deformation of the interface part which is inside
the defect not too close to its edges is $\eta=\delta p_{c}^{*}/\rho
g_{e}=\ell_{s}(\Pi/\Pi_{w}^{(1)})$, similar to that for the mesa defect
[namely Eq. (16) for $|x|<w/2$ except near the defect edges]. Thus, the
interface configuration before and after the jump is

$h_{\rm imb}^{-}=y_{d}$ (45)

and

$h_{\rm imb}^{+}=\begin{cases}y_{d},&|x|>w/2,\\\
y_{d}+\ell_{s}(\Pi/\Pi_{w}^{(1)}),&|x|<w/2.\end{cases}$ (46)

An exception to the validity of this calculation is for $\Pi/\Pi_{w}^{(1)}$
just above and very close to 1. In such case, $\Pi$ may be above
$\Pi_{w}^{(1)}$ yet below the exact threshold $\Pi_{c}^{(1)}$; even if not,
the influence of the edges extends very far inside the defect (see Fig. 3b).
For drainage (except for $\Pi/\Pi_{c}^{(1)}\approx 1$) immediately before the
jump the interface passes very close to the corners of the defect at
$y=y_{d}$; this is similar to the _mesa_ but not to the _regular_ narrow
defect case. It can be shown then that the interface shape around the defect
edge at $x=-w/2$ still has a center of symmetry, but at height
$y_{d}+\ell_{s}/2$, instead of in the corner of the defect. Mathematically,
$h(x_{0}+\Delta x)-y_{d}-\ell_{s}/2\approx-[h(x_{0}-\Delta
x)-y_{d}-\ell_{s}/2]$ when $|\Delta x|\ll w$, where $x_{0}$ is such that
$0<x_{0}+w/2\ll w$, and analogously for $x$ near $w/2$. The interfacial
configurations before and after the jump in drainage are then

$h_{\rm
dra}^{-}=\begin{cases}y_{d}-\ell_{s}{\displaystyle\frac{\Pi/\Pi_{w}^{(1)}-1}{2}},&|x|>w/2,\\\\[8.5359pt]
y_{d}+\ell_{s}{\displaystyle\frac{\Pi/\Pi_{w}^{(1)}+1}{2}},&|x|<w/2.\end{cases}$
(47)

and

$h_{\rm dra}^{+}=y_{d}-\ell_{s}\frac{\Pi/\Pi_{w}^{(1)}-1}{2}.$ (48)

These analytical expressions are confirmed by the numerical computations
(Appendix B), see Fig. 3b.

##  IV Energy balance: Dissipation during jumps

We begin by considering a closed imbibition-drainage cycle for a strong
defect. The first dissipation event occurs when the interface enters the
defect in imbibition. The external pressure is $\rho g_{e}h_{0}=\rho
g_{e}h_{\rm imb}^{-}=\rho g_{e}y_{d}$, and the work is $W_{\rm imb}=\rho
g_{e}y_{d}\int\limits_{-\infty}^{\infty}dx\,\eta_{\rm imb}(x)$; here we used
the fact that the deformation after the jump, $\eta_{\rm imb}(x)$, equals
$h_{\rm imb}^{+}-h_{\rm imb}^{-}$. The internal energy change during this
deformation is

$\displaystyle\Delta U_{\rm imb}$ $\displaystyle=U\\{h_{\rm
imb}^{+}\\}-U\\{h_{\rm imb}^{-}\\}$
$\displaystyle=\int\limits_{-\infty}^{\infty}dx\left[\frac{\gamma}{2}\left(\frac{d\eta_{\rm
imb}}{dx}\right)^{2}+\frac{\rho g_{e}}{2}(\eta_{\rm imb}^{2}+2y_{d}\eta_{\rm
imb})-\delta p_{c}^{*}\int_{y_{d}}^{y_{d}+\eta_{\rm imb}}dy\,F(x,y)\right],$
(49)

where $U\\{h\\}$ denotes the functional dependence of the internal energy on
the interface configuration. Thus, we obtain for the dissipated energy

$\Psi_{\rm imb}=\Delta U_{\rm imb}-W_{\rm
imb}=\int\limits_{-\infty}^{\infty}dx\left[\frac{\gamma}{2}\left(\frac{d\eta_{\rm
imb}}{dx}\right)^{2}+\frac{\rho g_{e}}{2}\eta_{\rm imb}^{2}-\delta
p_{c}^{*}\int_{y_{d}}^{y_{d}+\eta_{\rm imb}}dy\,F(x,y)\right].$ (50)

The second dissipation event occurs when the interface leaves the defect in
drainage. The work done on the interface is $W_{\rm dra}=-\rho g_{e}h_{\rm
dra}^{+}\int\limits_{-\infty}^{\infty}dx\,\eta_{\rm dra}(x)$, where the
deformation of the interface before the jump, $\eta_{\rm dra}(x)=h_{\rm
dra}^{-}-h_{\rm dra}^{+}$. The change of internal energy is

$\displaystyle\Delta U_{\rm dra}$ $\displaystyle=U\\{h_{\rm
dra}^{+}\\}-U\\{h_{\rm dra}^{-}\\}$
$\displaystyle=-\int\limits_{-\infty}^{\infty}dx\left[\frac{\gamma}{2}\left(\frac{d\eta_{\rm
dra}(x)}{dx}\right)^{2}+\frac{\rho g_{e}}{2}(\eta_{\rm dra}^{2}+2h_{\rm
dra}^{+}\eta_{\rm dra})-\delta p_{c}^{*}\int_{h_{\rm dra}^{+}}^{h_{\rm
dra}^{+}+\eta_{\rm dra}}dy\,F(x,y)\right].$ (51)

This provides the following dissipated energy:

$\displaystyle\Psi_{\rm dra}=\Delta U_{\rm dra}-W_{\rm
dra}=-\int_{-\infty}^{\infty}dx\left[\frac{\gamma}{2}\left(\frac{d\eta_{\rm
dra}(x)}{dx}\right)^{2}+\frac{\rho g_{e}}{2}\eta_{\rm dra}^{2}-\delta
p_{c}^{*}\int_{h_{\rm dra}^{+}}^{h_{\rm dra}^{+}+\eta_{\rm
dra}}dy\,F(x,y)\right].$ (52)

###  IV.1 Energy Dissipation: Mesa defect

For imbibition in a mesa defect, in the rightmost term in Eq. (50)

$\int_{y_{d}}^{y_{d}+\eta_{\rm imb}}dy\,F(x,y)={\rm Rect}(x/w)\eta_{\rm
imb}(x),$ (53)

which, together with the fact that $\eta_{\rm imb}$ equals $\eta$ from Eq.
(16), allows us to calculate the integral in Eq. (50), giving

$\Psi_{\rm imb}=-\frac{1}{2}\frac{\delta p_{c}^{*2}\ell_{c}}{\rho
g_{e}}\left[\frac{w}{\ell_{c}}+\exp(-w/\ell_{c})-1\right].$ (54)

Similarly, for drainage we substitute in Eq. (52) $h_{\rm dra}^{+}$ for
$h_{c}$ from Eq. (18), to obtain

$\int_{h_{\rm dra}^{+}}^{h_{\rm dra}^{+}+\eta_{\rm dra}}dy\,F(x,y)={\rm
Rect}(x/w)[\eta_{\rm dra}(x)+h_{c}-y_{d}].$ (55)

Using $\eta$ from Eq. (16) for $\eta_{\rm dra}$ in Eq. (52) gives

$\Psi_{\rm dra}=-\frac{1}{2}\frac{\delta p_{c}^{*2}\ell_{c}}{\rho
g_{e}}\left[1-\left(1+\frac{w}{\ell_{c}}\right)\exp(-w/\ell_{c})\right].$ (56)

The total dissipated energy for a closed hysteresis cycle, $\Psi_{\rm
tot}=\Psi_{\rm imb}+\Psi_{\rm dra}$, is given by

$\Psi_{\rm tot}=-\frac{w\delta p_{c}^{*2}}{2\rho
g_{e}}\left[1-\exp(-w/\ell_{c})\right].$ (57)

In the narrow defect limit, $w\ll\ell_{c}$, we get

$\Psi_{\rm imb}\approx\Psi_{\rm dra}\approx-\frac{1}{4}\frac{\delta
p_{c}^{*2}w^{2}}{\rho g_{e}\ell_{c}}.$ (58)

This dependence on system parameters is expected when analyzing the terms in
Eqs. (54) and (56) separately, as they are all of the same order of magnitude.
In particular, the dependence on the defect width $w$ is quadratic, as
$\eta_{m}\sim w$ [fixing all other parameters in Eq. (58)] and the width of
the region where $\eta\approx\eta_{m}$ is $w$-independent.

For the wide-defect limit, $w\gg\ell_{c}$, dissipation in imbibition reduces
to

$\Psi_{\rm imb}\approx-\frac{1}{2}\frac{\delta p_{c}^{*2}w}{\rho g_{e}}.$ (59)

Here, the dependence on $w$ is linear and matches those of the second and
third terms in Eq. (54); this is because $\eta_{\rm imb}$ in the wide-defect
limit is (i) approximately constant in the region of width $\approx w$, and
(ii) that constant is $w$-independent [fixing all other parameters in Eq.
(59)]. The first term in Eq. (54) is negligible (as the integrand peaks near
the edges of the defect in regions of width $\sim\ell_{c}\ll w$).

For drainage in the wide-defect limit, the dissipation is

$\Psi_{\rm dra}\approx-\frac{1}{2}\frac{\delta p_{c}^{*2}\ell_{c}}{\rho
g_{e}}.$ (60)

This is width-independent, and considerably smaller than in imbibition
($\Psi_{\rm imb}$), suggesting that the second and third terms approximately
cancel out. Thus, in the wide defect limit, dissipation can be thought of as
coming from the edges of the defect, rather than its whole width. This can
also be seen by examining the out-of-equilibrium pressure ($p_{e}$) during a
jump, which is used in Eqs. (8)–(9) to compute the dissipation as a force-
displacement product; this is illustrated in Video S1 in  Supplementary
Information (SI). The imbalance at the edges in drainage implies that it is
the source for a appreciable portion of the dissipation. We note that as the
example in Video S1 is for a moderately-wide defect (not the wide-defect
limit), a substantial part of the dissipation is associated with the motion of
the central parts of the interface.

###  IV.2 Energy Dissipation: Regular defect

For regular defects, the last term in the dissipation calculations, Eqs. (50)
and (52), becomes more complicated, because the defect consists of two parts.
With $F(x,y)$ given by Eq. (21), this term is

$\int_{y_{d}}^{y_{d}+\eta_{\rm imb}}dy\,F(x,y)={\rm
Rect}(x/w)\begin{cases}\eta_{\rm imb}^{2}/(2\ell_{s}),&\eta_{\rm
imb}<\ell_{s},\\\ \eta_{\rm imb}-\ell_{s}/2,&\eta_{\rm
imb}>\ell_{s},\end{cases}$ (61)

for imbibition, and

$\int_{h_{\rm dra}^{+}}^{h_{\rm dra}^{+}+\eta_{\rm dra}}dy\,F(x,y)={\rm
Rect}(x/w)\begin{cases}(\eta_{\rm dra}+h_{\rm
dra}^{+}-y_{d})^{2}/(2\ell_{s}),&\eta_{\rm dra}<\ell_{s}+y_{d}-h_{\rm
dra}^{+},\\\ \eta_{\rm dra}+h_{\rm dra}^{+}-y_{d}-\ell_{s}/2,&\eta_{\rm
dra}>\ell_{s}+y_{d}-h_{\rm dra}^{+}\end{cases}$ (62)

for drainage.

Calculating dissipation analytically (or even finding the expression for
$\eta$) for regular defects of an arbitrary width and slope is considerably
more difficult than for mesa defects. This is because the solution for the
interface can be in both the ramp and the plateau parts of the defect,
requiring matching between all the different parts of the interface. Numerical
results for $\eta_{\rm imb}$ and $\eta_{\rm dra}$ can be obtained as described
in Appendix B, followed by numerical integration to obtain $\Psi_{\rm imb}$
and $\Psi_{\rm dra}$. Nonetheless, analytical results can be obtained for
specific cases. First, as there is no dissipation for weak defects, the
dissipated energy for both imbibition and drainage vanishes as the defect
strength approaches the limit $\Pi/\Pi_{c}^{(1)}\to 1+$. On the other hand,
for $\Pi/\Pi_{c}^{(1)}\gg 1$, the interface solution crosses the defect
entirely within the plateau region, and therefore $\eta_{\rm imb}$ and
$\eta_{\rm dra}$ equal those for a mesa defect with the same $w$ and $\delta
p_{c}^{*}$. Therefore, the first two terms in Eqs. (50) and (52) are identical
for mesa and regular defects. If the interface jumps far into the plateau
region, the third term resembles that in the mesa case (e.g., for imbibition
$\eta_{\rm imb}-\ell_{s}/2\approx\eta_{\rm imb}$); if it is not (as is the
case for narrow defects in drainage, see Fig. 3a), the term is negligible.
Therefore, the dissipation for $\Pi/\Pi_{c}^{(1)}\gg 1$ is approximately equal
to that in a mesa defect (with the same $w$ and $\delta p_{c}^{*}$),

$\Psi_{\rm imb}\approx-\frac{1}{2}\frac{\delta p_{c}^{*2}\ell_{c}}{\rho
g_{e}}\psi_{\rm imb}=-\frac{\ell_{s}^{2}\ell_{c}^{3}\Pi^{2}}{2\gamma}\psi_{\rm
imb}$ (63)

and

$\Psi_{\rm
dra}\approx-\frac{\ell_{s}^{2}\ell_{c}^{3}\Pi^{2}}{2\gamma}\psi_{\rm dra}.$
(64)

Here,

$\psi_{\rm imb}=\frac{w}{\ell_{c}}+\exp(-w/\ell_{c})-1$ (65)

and

$\psi_{\rm dra}=1-\left(1+\frac{w}{\ell_{c}}\right)\exp(-w/\ell_{c}).$ (66)

Both the approach to zero dissipation as $\Pi\to\Pi_{c}^{(1)}$, and the
approach to the mesa results for $\Pi/\Pi_{c}^{(1)}\gg 1$, are confirmed by
Fig. 4, where numerical computations for a regular defect of an intermediate
width ($w/\ell_{c}=2$) are compared to analytical mesa results.

Figure 4: Dimensionless dissipated energy against dimensionless capillary
pressure perturbation for mesa and regular defects with $w/\ell_{c}=2$
($\ell_{s}/\ell_{c}=1$ for regular defects). For very strong defects (large
$\delta p_{c}^{*}$), the results for mesa and regular defects converge. The
difference between imbibition and drainage, which is significant for mesa
defects as well as for very strong regular defects, vanishes as $\delta p_{c}$
approaches the threshold between weak and strong.

Similarly, it is easy to calculate the dissipation for a regular defect with
an arbitrary slope $\Pi$, in the narrow defect limit ($w\ll\ell_{c}$). The
dissipation, that turns out to be equal in imbibition and drainage, is
obtained using Eqs. (37)–(40) for the interface shapes [$h_{\rm imb}^{-}$,
$h_{\rm imb}^{+}$, $h_{\rm dra}^{-}$ and $h_{\rm dra}^{+}$] in Eqs. (50),
(52), (61) and (62), as well as the expression for $\Pi_{n}^{(1)}$ given by
Eq. (34):

$\Psi_{\rm imb}=\Psi_{\rm
dra}=-\frac{w\ell_{s}^{2}}{2}\Pi\left(\frac{\Pi}{\Pi_{n}^{(1)}}-1\right)=-\frac{\ell_{s}^{2}\ell_{c}^{3}}{4\gamma}\Pi(\Pi-\Pi_{n}^{(1)})\left(\frac{w}{\ell_{c}}\right)^{2}.$
(67)

This equals the dissipation for the mesa defect in Eqs. (63)–(64) when both
(i) $w\ll\ell_{c}$, in which case we can approximate the expressions in Eqs.
(65)–(66) by $(1/2)(w/\ell_{c})^{2}$; and (ii) $\Pi\gg\Pi_{n}^{(1)}$, when
$\Pi(\Pi-\Pi_{n}^{(1)})\approx\Pi^{2}$ in Eq. (67).

It is useful to consider the following dimensionless quantity,

$\psi=-\frac{\Psi}{\ell_{s}^{2}\ell_{c}^{3}\Pi(\Pi-\Pi_{c}^{(1)})/(2\gamma)}.$
(68)

For arbitrary defect widths and slopes, Eqs. (27)–(29) predict that $\psi$
depends only on the dimensionless parameters $w/\ell_{c}$ and
$\Pi/\Pi_{c}^{(1)}$. According to Eq. (67), for narrow defects it is expected
to depend only on $w/\ell_{c}$ (and be identical for imbibition and drainage),
while for an arbitrary width $w$ it should approach the expressions in Eqs.
(65)–(66) as $\Pi/\Pi_{c}^{(1)}\to\infty$. This is demonstrated by plotting
the dissipation for various defect widths $w$ and $\Pi/\Pi_{c}^{(1)}$ (Fig.
5). Figure 5 also shows that as $\Pi/\Pi_{c}^{(1)}\to 1$, $\psi$ approaches a
finite value, and therefore the dissipation is linear in $\Pi-\Pi_{c}^{(1)}$
just above the threshold; this was shown in Eq. (67) for narrow defects, here
confirmed for an arbitrary width. This linear dependence is a consequence of
the finite interfacial jump; a faster approach to zero is expected when this
is not so. Moreover, the limits for $\psi$ as $\Pi/\Pi_{c}^{(1)}\to 1$ are the
same for imbibition and drainage, see also Fig. 4. Therefore, the dissipation
for imbibition and drainage is similar for narrow defects (regardless of
strength $\Pi$) and for arbitrary widths when
$\Pi-\Pi_{c}^{(1)}\ll\Pi_{c}^{(1)}$.

Figure 5:  Rescaled dissipated energy $\psi$ [Eq. (68)] in imbibition (blue)
and drainage (red) for various values of $\Pi/\Pi_{c}^{(1)}$ (a single regular
defect). For narrow defects ($w/\ell_{g}\ll 1$), the imbibition and drainage
curves collapse (independently of $\Pi/\Pi_{c}^{(1)}$), with $\psi\propto
w^{2}$. For arbitrary widths, for small slopes i.e. $\Pi/\Pi_{c}^{(1)}\to 1$,
imbibition and drainage approach the same limit, whereas as the slope
increases (“stronger” defect; black arrows), $\Pi/\Pi_{c}^{(1)}\to\infty$, the
dissipation approaches the bounding values of the dissipation in a mesa defect
(dashed gray), given by the expressions in Eqs. (65)–(66). For imbibition at
very wide defects, dissipation scales with the width, $\psi\propto w^{1}$.

The dependence of dissipation on the defect width for a regular defect is
similar to that for a mesa defect for fixed $\Pi\gg\Pi_{c}^{(1)}$. This can be
explained by considerations similar to those we used for the mesa defects. For
narrow defects, all the terms in Eqs. (50) and (52) are of the same order of
magnitude and $\propto w^{2}$. For wide defects, the second and third terms
are $\propto w$, giving rise to that dependence of $\Psi$ on $w$ in
imbibition, but they cancel out for drainage resulting in $w$-independent
dissipation. However, in contrast to mesa defects, for regular defects
dissipation also vanishes as $\Pi\to\Pi_{c}^{(1)}$, i.e. the terms that are of
the same order of magnitude must cancel out. This is easy to confirm for
imbibition in wide defects, where the first term is negligible and the last
two terms in the integrand are constant within the defect not too close to its
edges. The sum of these terms is

$\frac{\rho g_{e}}{2}\eta_{\rm imb}^{2}-\delta
p_{c}^{*}\int_{y_{d}}^{y_{d}+\eta_{\rm imb}}dy\,F(x,y)=\begin{cases}(1/2)(\rho
g_{e}-\Pi)\eta_{\rm imb}^{2},&\eta_{\rm imb}<\ell_{s},\\\ (1/2)(\rho
g_{e}-\Pi)\eta_{\rm imb}^{2}+(\Pi/2)(\eta_{\rm imb}-\ell_{s})^{2},&\eta_{\rm
imb}>\ell_{s},\end{cases}$ (69)

where we have used Eq. (61). This has the incorrect (positive) sign for any
$\eta_{\rm imb}>0$ when $\Pi<\Pi_{w}^{(1)}=\rho g_{e}$, thus, there are no
nonzero solutions for $\eta_{\rm imb}$, i.e., no stable deformed
configurations, and the defect is weak. For a strong defect
[$\Pi>\Pi_{w}^{(1)}$] and $\eta_{\rm imb}=\ell_{s}(\Pi/\Pi_{w}^{(1)})$ [see
Eqs. (45)–(46)], this sum is negative and approaches zero when
$\Pi\to\Pi_{w}^{(1)}$, as expected. We can also verify that the last two terms
in the expression for drainage cancel out for any $\Pi>\Pi_{w}^{(1)}$. For
$\eta>\ell_{s}+y_{d}-h_{\rm dra}^{+}$ the sum of these terms becomes

$\frac{\rho g_{e}}{2}\eta_{\rm imb}^{2}-\delta p_{c}^{*}\int_{h_{\rm
dra}^{+}}^{h_{\rm dra}^{+}+\eta_{\rm
dra}}dy\,F(x,y)=\frac{\Pi_{w}^{(1)}}{2}\eta_{\rm
dra}^{2}-\Pi\ell_{s}(\eta_{\rm dra}+h_{\rm dra}^{+}-y_{d}-\ell_{s}/2),$ (70)

where we used Eq. (62). Using Eq. (47) for $h_{\rm dra}^{-}=\eta_{\rm
dra}+h_{\rm dra}^{+}$ and Eq. (48) for $h_{\rm dra}^{+}$, we see that the sum
vanishes for all $\Pi$ values.

##  V cooperative origin of hysteresis and dissipation

A single mesa defect is strong, with interface displacement (e.g. PS
trajectory) which is hysteretic and irreversible (dissipative). A disordered
medium composed of multiple mesa defects, is also hysteretic and dissipative
[12]. From this, one might naively conclude that the overall hysteretic
response is simply the sum of the response of basic hysteretic entities (here,
“defects”). This conforms with the conventional thinking behind compartment
models such as Leverett and Preisach (where basic hysteretic entities are
called “hysterons”) [26]. Following the same logic, one would expect that a
medium composed of multiple weak defects—each of which is non-hysteretic,
reversible and non-dissipative—will also be non-hysteretic and reversible. To
examine this, we consider a simple system composed of two identical regular,
weak defects of width $w$ at a distance $d$ apart (measured between their
centers, with $d\geq w$), positioned at $y=y_{d}$ (Fig 1). Similarly to the
single defects considered in Sections. III and IV, each defect has a capillary
pressure profile which is linearly increasing in $y$ over a distance
$\ell_{s}$, after which it remains constant (plateau), with the total defect
length $\ell$. For each system with a given interdefect distance $d$, we
compute the energy dissipated in imbibition and drainage, analytically and
numerically, as described below.

###  V.1 Interactions between a pair of weak defects: Analytical evaluation

The local pressure balance in Eq. (10) is general, and thus it can be used for
any number of defects, regular or irregular. For the pair of weak defects
considered here, we write the capillary pressure profile as

$F(x,y)=a_{2}(x)c(y),$ (71)

where $c(y)$ has the same form as in the single regular defect case [Eq.
(23)], with the term $a_{1}(x)$ [Eq. (22)] replaced here with $a_{2}(x)$
consisting of two rectangular functions,

$a_{2}(x)={\rm Rect}\left(\frac{x+d/2}{w}\right)+{\rm
Rect}\left(\frac{x-d/2}{w}\right)=w\left[\delta_{w}(x+d/2)+\delta_{w}(x-d/2)\right],$
(72)

where we use $\delta_{w}(x)$ introduced in Section III.2.1. Note that for
$d=w$, the two-defect system reduces to a single defect of width $2w$.

Consider the case when the defects are narrow ($w/\ell_{c}\ll 1$), so that we
can replace $\delta_{w}(x)\to\delta(x)$. Using the resulting approximation of
$a_{2}(x)$ in (26) (substituting $a_{2}$ for $a_{1}$), under the assumption
that the interface solution crosses the defects within the ramps, the
deformation is

$\eta(x)=\Pi w[G_{0}(x+{d/2})+G_{0}(x-{d/2})][(h_{0}-y_{d})+\eta_{1}],$ (73)

where we set $\eta_{1}=\eta(-{d/2})=\eta({d/2})$ by using the symmetry of the
domain. Imposing $x={d/2}$ in (73), we obtain

$\eta_{1}=\frac{\Pi w\beta_{1}(h_{0}-y_{d})}{1-\Pi
w\beta_{1}}\quad\mbox{with}\quad\beta_{1}=G_{0}(d)+G_{0}(0)$ (74)

and therefore

$\eta(x)=\left[G_{0}(x+{d/2})+G_{0}(x-{d/2})\right]\frac{\Pi
w(h_{0}-y_{d})}{1-\Pi w\beta_{1}}.$ (75)

This is consistent with our assumption that the interface solution crosses the
defects within the ramps for sufficiently small positive $h_{0}-y_{d}$, when
the denominator is positive, i.e.,

$\Pi<\Pi_{n}^{(2)}=\frac{1}{w\beta_{1}}=\frac{1}{w\beta_{0}[\exp(-{d}/\ell_{c})+1]}=\frac{\Pi_{n}^{(1)}}{\exp(-d/\ell_{c})+1}\leq\Pi_{n}^{(1)}.$
(76)

As for the single defect, this is the criterion for weakness of the pair of
defects, and $\Pi_{n}^{(2)}$ is the (narrow-defect approximation of) the
critical slope for the transition from weak to strong. This method can be
extended to an arbitrary number of defects which can also have different
defect strengths (disordered media). Note that this critical slope for the
two-defect system is always smaller than the one for the single defect,
$\Pi_{n}^{(1)}$. This holds for defects of arbitrary width; for derivation of
the critical slope $\Pi_{c}^{(2)}$ for this case see Appendix A. For slopes
within $\Pi_{c}^{(2)}<\Pi<\Pi_{c}^{(1)}$, each defect is non-dissipative
(weak) when isolated (single defect only), whereas a system of two such
defects is dissipative and hysteretic (strong).

Within the narrow-defect approximation, let us consider now a situation where
the slope of the individual defects is $\Pi<\Pi_{n}^{(1)}$, that is, each
defect by itself is weak. The two-defect system is strong if
$\Pi>\Pi_{n}^{(2)}$, implying that

$\displaystyle\left[\exp(-d/\ell_{c})+1\right]>\frac{\Pi_{n}^{(1)}}{\Pi}.$
(77)

This inequality indicates that the two-defect system becomes strong if the
separation distance $d$ between the defect centers is $w\leq d<d_{c}$, where

$d_{c}=\ell_{c}\ln\left(\frac{\Pi}{\Pi_{n}^{(1)}-\Pi}\right).$ (78)

This is meaningful only for $d_{c}>w\approx 0$, that is, for
$\displaystyle\frac{\Pi_{n}^{(1)}}{2}\leq\Pi<\Pi_{n}^{(1)}$. Thus, there is a
factor of two between the largest and smallest slopes where a single defect is
weak but a pair can be strong. The factor of two reduces for wider defects; in
the limit $w/\ell_{c}\gg 1$, it approaches unity, as the critical slope
approaches $\rho g_{e}=\gamma/\ell_{c}^{2}$, independent of $d$.

Considering the case when the interface solution crosses the defects within
the plateaus, Eq. (25) (with $a_{1}$ replaced by $a_{2}$) gives for narrow
defects

$\eta(x)=\Pi w\ell_{s}[G_{0}(x+d/2)+G_{0}(x-d/2)].$ (79)

At the defects, this gives

$\eta(\pm d/2)=\frac{\Pi
w\ell_{s}\ell_{c}}{2\gamma}[1+\exp(-d/\ell_{c})]=\frac{\Pi}{\Pi_{n}^{(2)}}\ell_{s},$
(80)

a result analogous to the single-defect case [see Eq. 36]. For dissipation
calculations, since the width of the defect cannot be neglected, we write

$\eta(x)=\frac{\Pi}{\Pi_{n}^{(2)}}\ell_{s}\begin{cases}\exp\\{-[|x|-(d+w)/2]/\ell_{c}\\},&|x|>(d+w)/2,\\\
\cosh[x/\ell_{c}]/\cosh[(d-w)/(2\ell_{c})],&|x|<(d-w)/2,\\\
1,&(d-w)/2<|x|<(d+w)/2,\end{cases}$ (81)

where we have assumed that $\eta(x)$ is constant within the defects, and equal
to Eq. (80), and the first two lines are, essentially, Eq. (79), with $d$
replaced by $d+w$ in the first line and $d-w$ in the second (a negligible
change) to make the result continuous. If the pair is strong
($\Pi/\Pi_{n}^{(2)}>1$), during imbibition the interface experiences a jump
when it first touches the defects, just as in the single-defect case; thus,
the interface configurations before and after the jump are

$h_{\rm imb}^{-}=y_{d}$ (82)

and

$h_{\rm
imb}^{+}=y_{d}+\frac{\Pi}{\Pi_{n}^{(2)}}\ell_{s}\begin{cases}\exp\\{-[|x|-(d+w)/2]/\ell_{c}\\},&|x|>(d+w)/2,\\\
\cosh[x/\ell_{c}]/\cosh[(d-w)/(2\ell_{c})],&|x|<(d-w)/2,\\\
1,&(d-w)/2<|x|<(d+w)/2.\end{cases}$ (83)

Likewise, during drainage the interface behaves similarly to the single-defect
case, where the jump occurs when it reaches the boundary between the plateau
and the ramp, with

$h_{\rm
dra}^{-}=\begin{cases}y_{d}-\ell_{s}\left\\{\frac{\Pi}{\Pi_{n}^{(2)}}(1-\exp\\{-[|x|-\frac{d+w}{2}]/\ell_{c}\\})-1\right\\},&|x|>(d+w)/2,\\\
y_{d}-\ell_{s}\left\\{\frac{\Pi}{\Pi_{n}^{(2)}}(1-\cosh[\frac{x}{\ell_{c}}]/\cosh[\frac{d-w}{2\ell_{c}}])-1\right\\},&|x|<(d-w)/2,\\\
y_{d}+\ell_{s},&(d-w)/2<|x|<(d+w)/2\end{cases}$ (84)

and

$h_{\rm dra}^{+}=y_{d}-\ell_{s}(\Pi/\Pi_{n}^{(2)}-1).$ (85)

The corresponding calculations of interface shapes for defects of an arbitrary
width can be done numerically using the method in Appendix B.

The dissipated energy in imbibition and drainage can be evaluated using eqs.
(50) and (52), with the values of $h_{\rm imb}^{\pm}$ and $h_{\rm dra}^{\pm}$
computed analytically or numerically, and the expressions for the integrals of
$F(x,y)$ from Eqs. (61)–(62), where ${\rm Rect}(x/w)$ are replaced by ${\rm
Rect}([x+d/2]/w)+{\rm Rect}([x-d/2]/w)$. For the narrow defect approximation,
using Eqs. (82)–(85), as well as Eq. (76) for $\Pi_{n}^{(2)}$, the result
reads

$\Psi_{\rm imb}=\Psi_{\rm
dra}=-w\ell_{s}^{2}\Pi\left[\frac{\Pi}{\Pi_{n}^{(2)}}-1\right].$ (86)

This expression resembles its counterpart for the single-defect case, Eq.
(67), except for the critical threshold at which dissipation approaches zero
[$\Pi_{n}^{(2)}$ instead of $\Pi_{n}^{(1)}$] and the factor of two; this is
intuitive, as for two defects far apart the dissipation is additive, and the
threshold remains the same. Equation (86) is also consistent with the fact
that two touching defects ($d=w$) are equivalent to a single defect of twice
the width. Finally, we find that the accuracy of the narrow-defect theory can
be improved upon replacing in Eq. (86) the approximate threshold,
$\Pi_{n}^{(2)}$, with the exact result, $\Pi_{c}^{(2)}$ [Eqs. (90)–(91)],
providing

$\Psi_{\rm imb}=\Psi_{\rm
dra}=-wl^{2}\Pi\left[\frac{\Pi}{\Pi_{c}^{(2)}}-1\right].$ (87)

###  V.2 Interactions between a pair of weak defects: Numerical verification

To examine these intriguing theoretical predictions, we use numerical
computations and simulations varying the distance $d$ for a fixed defect shape
(slope $\Pi$, ramp length $\ell_{s}$, width $w$), computing the energy
dissipated during the imbibition and drainage trajectories (see Appendix C for
parameter values). We calculate the dissipated energy using (i) the numerical
computations described in Appendix B; (ii) direct numerical simulations of the
interface evolution and the corresponding energy dissipation [12] (referred to
as simulations, to distinguish from the numerical computations of Appendix B).
While the numerical simulations are more computational costly than the
computations in Appendix B, the simulations can be used for any arbitrary
capillary pressure field $p_{c}(x,y)$ (e.g. disorder with prescribed defect
strength distributions in [12]).

Our numerical evaluations show that for pairs of weak defects sufficiently far
apart, $d>d_{c}$ [where $d_{c}$ is approximated by Eq. (78)], there is no
dissipation and hysteresis (Video S2 in  SI), whereas for $d<d_{c}$,
dissipation emerges (Fig. 6). This dissipation arises from abrupt jumps of the
interface along the defect slope in both imbibition and drainage, such that
the PS response becomes hysteretic; the closer the defects are, the stronger
the dissipation (and the width of the hysteresis cycle, e.g. see Videos S3–S4
in  SI).

Figure 6:  Emergence of energy dissipation (“ _strong_ ” behavior) in a pair
of _weak_ defects. As long as the defects are sufficiently far apart
($d>d_{c}$), the system is reversible and there is no dissipation. Once
$d<d_{c}$, the displacement becomes irreversible, and dissipation and
hysteresis emerge. We plot the dissipation calculated from simulations ([12])
and numerical computation (Appendix B; dotted lines), where blue and red refer
to imbibition and drainage, respectively. Also plotted are the analytical
solutions from two variants of the narrow defect theory: (1) Eq. (86)
(approximate $\Pi_{n}^{(2)}$ [dash-dot]; and (2) the more accurate Eq. (87)
with an exact value of $\Pi_{c}^{(2)}$ [thick dashed].

As expected for narrow defects, [Eq. (86)], the numerical computations give
nearly identical dissipated amounts for imbibition and drainage. The
simulations provide similar results, though slightly exaggerate the difference
between imbibition and drainage, due to numerical discretization errors (can
be reduced by refinement, at the expense of higher computational cost). The
distance at which dissipation vanishes in the numerical computations is
identical to the exact result for the critical distance from Eqs. (90)–(91)
(“theory(2)” in legend of Fig. 6). The narrow defect theory [dissipation from
Eq. (86) and critical distance from Eq. (78); “theory(1)” in legend of Fig. 6]
provides a reasonable approximation. For wider defects, the deviation between
the dissipation computed in imbibition and drainage, and between these and the
dissipation evaluated with Eq. (87), increases. We note that the close
agreement between the simulated dissipation for imbibition and Eq. (86) is
coincidental.

###  V.3 Interactions between a pair of weak defects: Experimental observation

Here, we provide an experimental proof-of-concept showing the emergence of
dissipation in a pair of weak defects as they are brought close to each other.
We use 3-D printing to manufacture a series of systems (imperfect Hele-Shaw
cells), with a single weak defect and two pairs of identical defects at two
different separation distances. As we cannot measure the energy dissipated
experimentally, we use the hysteresis cycle as a proxy for reversibility; the
larger the area within a closed PS cycle is, the larger the dissipation [12].
Details of the experiments including the manufacturing, setup and image
analysis, are provided in Appendix D.

Our experiments validate the findings we obtained theoretically and
numerically: while a single regular defect of given geometry (weak) can behave
reversibly, showing no hysteresis (and thus no dissipation) (Fig. 7a; see also
Video S5 in  SI), a pair of defects (each of identical geometry to the former)
close enough together becomes hysteretic (Fig. 7b; Video S6 in  SI) due to the
spatial interactions between the otherwise reversible entities. Decreasing the
pair separation increases the hysteresis (Fig. 7c; Video S7 in  SI). In Fig.
7, we measure the maximum deformation $\eta_{m}$ along the middle line of the
defect, and the baseline position $h_{f}$ as the vertical distance between the
unperturbed interface (far from the defect) and the bottom of the defect,
i.e., $h_{f}=h_{0}-y_{d}$ (see also Fig. 1).

Figure 7:  Imbibition-drainage cycle for a single regular defect (a) and a
pair of regular defects of width $w=10$ mm, with a separation (measured
between their centres) $d$ of (b) 13 mm and (c) 11 mm, respectively. For each
case, we show the imbibition-drainage trajectories in terms of the interface
deformation $\eta_{m}$ vs. its baseline position $h_{f}$ (left panels) and an
experimental image of the interface (right). Here $\eta_{m}$ is the distance
between the blue dot and the red line. Defects are highlighted by green
rectangles. While the interface moves through a single regular defect
reversibly, i.e. with no hysteresis and therefore no dissipation (a), a pair
of such defects close enough together exhibits hysteresis (b), which increases
as their separation distance decreases (c), confirming our theoretical
prediction.

We note that in these experiments, effects that are not considered in our 2-D
linear model can be of importance 111 We also note that our 2-D model is
strictly valid only when the out-of-plane thickness varies slowly in space,
while in the examples considered here the defects contain abrupt changes in
thickness.. One such effect is large interface deformation resulting in
nonlinear curvature, where in our 2-D model we (i) approximate the total
curvature as the sum of in- and out-of-plane components, and (ii) use a linear
approximation for the former (vs. the full nonlinear formulation, e.g. see
[23]). Further differences between our 2-D model and the experiments arise
from 3-D effects related to (i) the curvature of the meniscus between adjacent
defects; (ii) the finite width of the meniscus in-plane projection; and (iii)
thin liquid film advancing ahead of the experimentally-observed _2-D
projection_ of the (3-D) meniscus, resulting in an uncertainty in its
position, hence $\eta_{m}>0$ at $h_{f}<0$ in Fig. 7; see further elaboration
in Appendix D. These effects preclude a quantitative comparison between our
theoretical predictions and the experiments, which are presented here as
proof-of-concept that qualitatively supports our theory.

##  VI Conclusions

In this paper, we study energy dissipation during quasi-static fluid-fluid
displacements across an imperfect Hele-Shaw cell (with “defects” i.e. local
perturbations in thickness and thus capillary pressure $p_{c}$). To explore
the fundamental mechanisms for energy dissipation and hysteresis, we consider
a simple system comprising of isolated (single) defects of two shapes:
“regular” with _continuous_ (here linear) variations in $p_{c}$ (in the
direction of fluid advancement), and “mesa” defects where $p_{c}$ changes
_abruptly_. Depending on the slope of $p_{c}$, the system can be classified as
either (i) “weak”, where the interface passage across the defect (imbibition
and drainage) is fully reversible with no dissipation and hysteresis; and (ii)
“strong”, with dissipation and hysteresis. We derive analytical expressions
for the evolution of the fluid-fluid configuration as it deforms when it
passes through defects and the corresponding energy dissipation, which are
validated numerically. The analysis exposes subtle morphological differences
between imbibition and drainage, originating from the differences in wet area
when a jump occurs.

A novel intriguing finding is that irreversible, hysteretic behavior and the
associated energy dissipation can emerge from the interaction of weak
defects—objects which are intrinsically (when isolated) non-dissipative and
non-hysteretic. This is demonstrated for a simple system: a pair of identical
weak defects, varying their separation distance. While far enough apart, the
pair of defects behaves as weak. Once the defect distance falls below a
threshold—dictated by surface tension and the system parameters (defect width
and slope of $p_{c}$), it becomes strong, producing hysteresis and
dissipation. We compute this threshold and the energy dissipated analytically,
for the approximation of narrow defects (width much smaller than the capillary
length), and numerically (for arbitrary width) using two methods: (i)
numerical simulations of the sequence of equilibrium configurations and (ii)
numerical computation of the equilibrium conditions. A proof-of-concept
experiment using 3-D printed cells demonstrates the emergence of irreversible,
hysteretic behavior due to the interaction between non-hysteretic, reversible
defects, validating our theoretical findings.

The cooperative mechanism exposed here is uniquely highlighted by our model.
In contrast, classical models for hysteresis such as compartment models [26]
rely on the existence of a basic hysteretic unit (“hysteron”), and the overall
hysteretic response is simply the sum of the responses of these units. In the
zero-temperature RFIM isolated spin flips are reversible [16], and thus
individual spins play a role equivalent to our weak defects. Similarly to the
present case, hysteresis and dissipation in the RFIM emerge from the
collective response of the system, in the form of spin-flip avalanches
triggered by the local spin-spin interactions. A crucial difference, however,
is that spin-spin interactions in RFIM are present in the whole system, and
thus irreversible events can occur anywhere. In contrast, the interaction
between defects in our system is mediated by the two-phase interface, so that
irreversible jumps and dissipation are localized at the interface.
Interestingly, the emergence of hysteresis as a collective response of a
system of individually non-hysteretic agents with continuum responses,
interacting at each step through the memory of the predecessors, was applied
to explain irreversibility in climate events [32].

To the best of our knowledge, our model is the first to describe how
hysteresis of pressure-saturation curves during cyclic flows in disordered
media emerges from the interactions among defects, in a physical, non-
phenomenological manner. Specifically in the context of fluid displacement,
both (i) the “ink bottle effect”—the canonical conceptual model for the pore-
scale mechanism underlying pressure-saturation hysteresis [3], as well as (ii)
the “Leverett model”—a compartment model predicting macroscopic hysteresis
[26], do not include the key concept of spatial distance among the basic units
(pores or defects). An oversimplification of the interactions between basic
model units is also inherent to the classical capillary bundle model used in
soil physics to predict the pressure-saturation relationship [33]. We argue
that the cooperative mechanism unveiled here is the dominant mechanism
responsible for hysteresis in multiphase fluid systems, noting that other
mechanisms such as contact angle hysteresis in surface wetting [34, 35],
changes in interfacial connectivity (snap-off) and fluid trapping [2] also
contribute to hysteresis.

In conclusion, we present a detailed, rigorous investigation of the
fundamental process of energy dissipation between consecutive metastable
configurations in the passage of an interface through topographical defects.
In the quasi-static limit considered here, viscous dissipation due to finite
velocity of the interface displacement is disregarded. The simplicity of this
conceptual model allows systematic theoretical examination of the origins of
energy dissipation and hysteresis. The insights gained here considering
individual defects are of immediate relevance to disordered media containing
multiple interacting defects [12]. Interesting directions for further studies
are the introduction of mechanisms that are not considered in the current
model system, to gain understanding of two-phase flow in more complex
disordered systems such as porous and fractured materials [8]. One is dynamic
effects associated with viscous dissipation in rate-driven systems and at high
flow rates, connecting insights gained from detailed investigation of single
isolated defects [23]. Another is the nonlinear, 3-D effects that were evident
even in the simple experiments presented here. Finally, links between our
novel modelling approach and the RFIM and other lattice models open the way to
gain fundamental understanding of how cooperative interactions among non-
hysteretic, non-dissipative entities could lead to the emergence of hysteresis
and dissipation in a wide variety of driven disordered systems [36, 37, 38].

##  Appendix A Evaluating the critical slope

Here, we derive the critical pressure slope for a _pair_ of identical defects,
from which we obtain the result for a single defect by putting the two defects
next to each other (equivalent to a single defect of twice the width), or
infinitely far apart (equivalent to a single defect of the same width). To
analyze defects of an arbitrary width, we use the original differential
equation for the interface shape, Eq. (1), instead of the Green’s function-
based approach.

Consider a pair of defects of width $w$, with distance $d$ between their
symmetry axes, as introduced in Section V. The capillary pressure profile is
given by Eqs. (71), (72) and (23). For the interface profile with
$h_{0}=y_{d}$, if the interface solution crosses the defects entirely within
the ramp regions, the equation for $\eta(x)=h(x)-y_{d}$ is

$\gamma\eta^{\prime\prime}=\begin{cases}(\rho g_{e}-\Pi)\eta,&||x|-d/2|<w/2\
\text{(inside the defects)},\\\ \rho g_{e}\eta,&||x|-d/2|>w/2\ \text{(outside
the defects)}.\end{cases}$ (88)

For $\Pi>\rho g_{e}$ the solutions are given by

$\eta(x)=\begin{cases}C_{1}\cosh(x/\ell_{c}),&|x|<(d-w)/2,\\\
C_{2}\cos[\sqrt{\Pi/\gamma-1/\ell_{c}^{2}}x]+C_{3}\sin[\sqrt{\Pi/\gamma-1/\ell_{c}^{2}}|x|],&(d-w)/2<|x|<(d+w)/2,\\\
C_{4}\exp(-|x|/\ell_{c}),&|x|>(d+w)/2.\end{cases}$ (89)

The top line in (89) corresponds to the space between the defects, the middle
to the part of the interface inside the defects, and the bottom to the outer
regions. Matching conditions, requiring that the values of $\eta(x)$ and its
first derivative $\eta^{\prime}(x)$ match at $|x|=(d-w)/2$ and $|x|=(d+w)/2$,
provide a system of four linear homogeneous equations for four unknowns
$C_{1}$–$C_{4}$, which has nontrivial solutions when its determinant is zero.
Given the parameters $w$, $d$, $\ell_{c}$ and $\gamma$, this can only be
satisfied for a single value of $\Pi=\Pi_{c}^{(2)}$. For this value, a
continuum of nontrivial solutions exist, differing by the value of the
prefactor from zero to the maximum value for which the solution still crosses
the defects entirely within the ramps.

The interpretation of this result is as follows. (i) For $\Pi<\Pi_{c}^{(2)}$,
only the trivial solution $\eta(x)=0$ of the original equation (1) with
$h_{0}=y_{d}$ exists; as $h_{0}$ increases, this solution evolves
continuously, thus, the defect is weak. (ii) As $\Pi$ increases, the rate of
deformation, $d\eta_{m}/dh_{0}$ grows, approaching infinity as
$\Pi\to\Pi_{c}^{(2)}$ from below, so that at $\Pi=\Pi_{c}^{(2)}$ an
infinitesimal change in $h_{0}$ gives rise to a finite deformation. This is
consistent with the existence of a continuum of solutions with $h_{0}=y_{d}$
at $\Pi=\Pi_{c}^{(2)}$. (iii) for $\Pi>\Pi_{c}^{(2)}$ there are no nontrivial
solutions that cross the defects within the ramps, but there is one at least
partially within the plateaus and thus a jump occurs at $h_{0}=y_{d}$ and the
defect is strong. Then $\Pi_{c}^{(2)}$ is the critical value separating weak
and strong defect pairs.

By equating the determinant of the above-mentioned system of four equations to
zero, an equation for $\Pi_{c}^{(2)}$ can be obtained. We define a quantity
$\zeta$ such that

$\Pi_{c}^{(2)}=\gamma[1/\ell_{c}^{2}+(\zeta/w)^{2}].$ (90)

Then, the following equation

$\tan\zeta=\frac{1+\tanh[(d-w)/(2\ell_{c})]}{(\ell_{c}\zeta/w)^{2}-\tanh[(d-w)/(2\ell_{c})]}\frac{\ell_{c}\zeta}{w}$
(91)

has one solution in the interval $0<\zeta<\pi$, which, generally, needs to be
found numerically, and then $\Pi_{c}^{(2)}$ is given by Eq. (90). If $\Pi$,
$w$, $\ell_{c}$ and $\gamma$ are given, a fully analytical solution for the
critical value of the distance $d$ is possible.

Next, we verify that for narrow defects ($w/\ell_{c}\ll 1$), the result of Eq.
(76) is recovered. Assuming (to be confirmed by the calculation) that
$w/\ell_{c}\ll\zeta\ll 1$, we can approximate Eq. (91) as

$\zeta=\\{1+\tanh[d/(2\ell_{c})]\\}\frac{w}{\ell_{c}\zeta},$ (92)

the solution of which is

$\zeta=\left[\frac{w}{\ell_{c}}\left(1+\tanh\frac{d}{2\ell_{c}}\right)\right]^{1/2}.$
(93)

This indeed satisfies the above inequality for $\zeta$. Then, we write Eq.
(90),

$\Pi_{c}^{(2)}=\gamma\left[\frac{1}{\ell_{c}^{2}}+\frac{1+\tanh(d/2\ell_{c})}{w\ell_{c}}\right]\approx\frac{\gamma[1+\tanh(d/2\ell_{c})]}{w\ell_{c}},$
(94)

which, after a simple transformation, coincides with Eq. (76). On the other
hand, for $w\gg\ell_{c}$, since $\zeta$ is finite, Eq. (90) gives
$\Pi_{c}^{(2)}=\gamma/\ell_{c}^{2}$ for any $d$.

Finally, we obtain the critical slope for a single defect by noting that for
$d=w$ (two defects put together with no gap), we get

$\tan\zeta=\frac{w}{\ell_{c}\zeta}.$ (95)

This corresponds to a single defect of width $2w$; then, for a single defect
of width $w$,

$\tan u=\frac{w}{2\ell_{c}u},$ (96)

and

$\Pi_{c}^{(1)}=\gamma[1/\ell_{c}^{2}+(2u/w)^{2}].$ (97)

Equation (96) is a transcendental equation for $u$ that needs to be solved
numerically. However, a fully analytic solution is possible for the critical
value of $w$ given $\Pi$, $\ell_{c}$ and $\gamma$. We can also check that for
two defects very far apart ($d-w\gg\ell_{c}$) this single-defect result is
recovered. Indeed, in this case Eq. (91) becomes

$\tan\zeta=\frac{2\ell_{c}\zeta/w}{(\ell_{c}\zeta/w)^{2}-1}=\frac{2(w/\ell_{c}\zeta)}{1-(w/\ell_{c}\zeta)^{2}}.$
(98)

Using the trigonometric identity $\tan 2X=2\tan{X}/(1-\tan^{2}X)$, we get

$\tan(\zeta/2)=\frac{w}{\ell_{c}\zeta},$ (99)

which coincides with Eq. (96) if $u=\zeta/2$; Eq. (97) then coincides with Eq.
(90).

##  Appendix B Mixed numerical-analytical computation of the interface profile

Here, we describe the mixed numerical-analytical method we have used to
calculate the interface configuration $h(x)$, in particular, (i) after the
jump during imbibition and (ii) before the jump during drainage. As similar
approaches have been used for a single and a pair of defects, we describe both
cases at the same time, indicating differences where applicable.

We consider capillary pressure profiles given by Eqs. (21)–(23) for a single
defect, and Eqs. (71), (72), (23) for a pair. Thus, the defect consists of a
ramp (slope) of length $\ell_{s}$, followed by a plateau. We only consider
cases where in the range(s) of $x$ where the defect(s) is (are) located, the
interface is entirely within the defect(s), i.e. for $|x|<w/2$ (single defect)
or $||x|-d/2|<w/2$ (pair), $y_{d}<h(x)<y_{d}+\ell$. The interface then obeys
the following equation for a single defect,

$\gamma h^{\prime\prime}=\begin{cases}\rho g_{e}(h-h_{0}),&|x|>w/2,\\\ \rho
g_{e}(h-h_{0})-\Pi(h{-y_{d}}),&|x|<w/2\ \text{and}\ h<y_{d}+\ell_{s},\\\ \rho
g_{e}(h-h_{0})-\Pi(\ell_{s}{-y_{d}}),&|x|<w/2\ \text{and}\
h>y_{d}+\ell_{s}\end{cases}$ (100)

and for a defect pair,

$\gamma h^{\prime\prime}=\begin{cases}\rho g_{e}(h-h_{0}),&|x|>(d+w)/2,\
|x|<(d-w)/2\\\ \rho g_{e}(h-h_{0})-\Pi(h{-y_{d}}),&||x|-d/2|<w/2\ \text{and}\
h<y_{d}+\ell_{s},\\\ \rho g_{e}(h-h_{0})-\Pi(\ell_{s}{-y_{d}}),&||x|-d/2|<w/2\
\text{and}\ h>y_{d}+\ell_{s}.\end{cases}$ (101)

In the outer region ($|x|>w/2$ for a single defect and $|x|>(d+w)/2$ for a
pair), the solution is

$h=h_{0}+\eta_{b}\exp(-\Delta x/\ell_{c}),$ (102)

where $\Delta x=|x|-w/2$ for a single defect or $\Delta x=|x|-(d+w)/2$ for a
pair. The constant $\eta_{b}$, the value of $h$ at the (outer) boundary of the
defect, needs to be found based on the requirement that the solution is
symmetric, thus, $h^{\prime}(0)=0$, by matching to other parts of the
solution, as discussed below; the solution for $\eta_{b}$ may or may not exist
depending on $h_{0}$ and the parameters of the defect(s).

Inside the defect(s) ($|x|<w/2$ for a single defect or $||x|-d/2|<w/2$ for a
pair), the solution, in general, consists of pieces of functions that can be
found analytically. As the matching between these pieces is cumbersome, that
part of the solution is obtained here by numerical integration. Assuming that
$h_{0}$ and $\eta_{b}$ are known, the values of $h$ and its first derivative
are provided from the outer solution at the $x=w/2$ [or $x=(d+w)/2$] boundary
of the defect, and can serve as the initial conditions for numerical
integration. For a single defect, integration can be carried out down to $x=0$
to find out if the $h^{\prime}(0)=0$ condition is satisfied. Thus, we can find
$\eta_{b}$ by solving the $h^{\prime}(0)=0$ equation using the bisection root-
finding scheme [39, Chapter 9.1]. For a pair of defects, $x=0$ is in the
middle between the defects; in that region between the defects, the analytical
solution with $h^{\prime}(0)=0$ is

$h(x)=h_{0}+C\cosh(x/\ell_{c}),$ (103)

where $C$ is an unknown constant. This gives the condition on the values of
$h$ and its first derivative on the inner boundary of the defect,

$\frac{h^{\prime}([d-w]/2)}{h([d-w]/2)-h_{0}}=\frac{1}{\ell_{c}}\tanh\frac{d-w}{2\ell_{c}}.$
(104)

By integrating numerically down to $(d-w)/2$ within the defect, we find if
this condition is satisfied, which, as before provides an equation for the
bisection scheme to find $\eta_{b}$. We note that for $\Pi>\rho g_{e}$ the
solution for $|x|<w/2$ (or $||x|-d/2|<w/2$) can be oscillatory and it is
possible that several roots $\eta_{b}$ and associated solutions exist.
However, apart from the solution with a single maximum inside (each) defect,
these solutions go below $y_{d}$ and as such violate the conditions outlined
above.

During imbibition, the interface jump happens once the interface touches the
defect. Thus, $h_{\rm imb}^{-}=y_{d}$ and $h_{\rm imb}^{+}$ is the solution of
Eq. (100) or (101) with $h_{0}=y_{d}$. This solution can be found as described
above. For drainage, the situation is more complicated, because $h_{0}$ is
unknown. Its value is a bifurcation point such that for lower values there are
no solutions crossing the defect. In other words, for lower $h_{0}$,
regardless of the value of $\eta_{b}$, the defect will not sufficiently deform
the interface. For a single defect this means that the maximum possible value
of $h^{\prime}(0)$ is below zero; for a pair, the maximum possible value of
the ratio on the left-hand side of Eq. (104) is smaller than its right-hand
side. We then carry out a nested procedure, where in the inner cycle, for a
particular $h_{0}$ we find the maximum value of $h^{\prime}(0)$ (for a single
defect) or of the left-hand side of Eq. (104) (for a pair), using the golden-
section algorithm [39, Chapter 10.2], and then, via bisection, find the value
of $h_{0}$ for which this maximum value is zero (for a single defect) or the
right-hand side of Eq. (104) (for a pair). This eventually provides both
$h_{\rm dra}^{-}$ and $h_{\rm dra}^{+}$ (the latter equal to the value of
$h_{0}$ resulting from the procedure).

The above procedure was used to find the interface shapes for narrow and wide
single defects in Fig. 3. For clarity, we add a similar plot for the general
case of a defect of an intermediate width ($w/\ell_{c}=2$; Fig. 8). Overall,
qualitatively, the results are intermediate between the limits of wide and
narrow defects. However, notably, for imbibition, while in both limits
$\eta(0)/\ell_{s}=\Pi/\Pi_{c}^{(1)}$, here the value of $\eta(0)$ is slightly
higher, thus, the dependence of $\eta(0)$ on the defect width is non-
monotonic.

Figure 8:  Interface profiles before and after jumps for an intermediate-width
($w/\ell_{c}=2$) regular defect, for different capillary pressure slopes. The
ramp part of the defect is shown as gradient shading and the plateau part in
uniform gray. The boundary between these two regions at $h/\ell_{c}=1$ is
marked with a green dotted line. Dashed and solid lines represent the
undeformed (flat) and deformed configurations, i.e. before and after the jump
in imbibition, and the opposite in drainage. These results are complementary
to those shown in Fig. 3.

Once the interface profiles before and after the jump are found, the
dissipated energy can be calculated using Eqs. (50) and (52), where the
integration is carried out numerically inside the defect(s) and analytically
outside. Note that if one fixes $\eta_{b}$ (for imbibition) or $h_{0}$ (for
drainage), the corresponding value of $d$ for which Eq. (104) is satisfied can
be found directly, avoiding bisection. This can be used to speed up
computations when obtaining $\Psi$ vs $d$ data (Fig. 6), by varying $\eta_{b}$
(or $h_{0}$) and producing a table of [$\Psi$, $d$] values.

##  Appendix C Numerical details

For the results shown in Fig. 6 (and Videos S2–S4 in  SI), we use the
following parameters: $g_{e}=0.2$ m/s2, $\rho=998$ kg/m3, $\gamma=20.7$ mN/m
(such that $\ell_{c}\approx 10.2$ mm). The defect width is $w=0.8$ mm, which
means that it is relatively narrow ($w/\ell_{c}\approx 0.08$), and we expect
the narrow defect theory to yield a good approximation. The defect profile is
such that $\delta p_{c}$ changes linearly between zero and $\approx$$7.39$ Pa
along a slope of length $\ell_{s}=2$ mm. For these values, an isolated single
defect is weak, as the pressure slope, $\Pi\approx 3.70\times 10^{3}$ Pa/m, is
smaller than the critical value for a single defect, computed using both (i)
the narrow-defect approximation, with Eq. (34) giving $\Pi_{n}^{(1)}\approx
5.08\times 10^{3}$ Pa/m; and (ii) the exact arbitrary-width result of Appendix
A, with Eqs. (96)–(97), such that $\Pi_{c}^{(1)}\approx 5.22\times 10^{3}$
Pa/m.

##  Appendix D Experimental details

An imperfect Hele-Shaw cell was manufactured using stereolithography 3-D
printing. The cell was produced with three sets of defects: a single defect
and two sets of defect pairs with varying inter-defect gaps, see Fig. 9 (a).
This efficient design allows us to run three separate experiments on the same
cell in sequence.

Figure 9: (a) A 3-D printed including the 3 sets of defects (left), with a
typical image of the interface during an experiment (right), showing the
interface deformation due to its interaction with the top defects (dark green
rectangles). Silicone oil enters from the bottom to displace the air initially
filling the cell. The oil-air interface advances (imbibition) and recedes
(drainage) in response to an increase or decrease in the oil pressure at the
bottom of the cell, respectively. The cell is tilted by an angle
$\alpha=38^{\circ}$ with respect to the horizontal. (b) Imbibition-drainage
cycle for the single defect. For each frame analyzed, we track the unperturbed
position of the interface far from the defect $h_{f}$ and the maximum
perturbation $\eta_{m}$ of the interface with regards to the unperturbed
position. The presented cycle is reversible, namely the defect is _weak_. In
panel (b), the gray shaded area corresponds to the data shown in Fig. 7 (a),
where the arrows point to interesting physical effects not considered in our
theoretical analysis (see text).

To experimentally produce variations in the capillary pressure $\delta p_{c}$,
we introduce modulations of the Hele-Shaw cell gap space $b(x,y)=b_{0}-\delta
b(x,y)$. This introduces a local variation of the out-of-plane capillary
pressure $p_{c}(x,y)=p_{c}^{0}+\delta p_{c}(x,y)$, where
$p_{c}^{0}=2\gamma\cos(\theta)/b_{0}$ ($\theta$ being the contact angle) and
$\delta p_{c}=p_{c}^{0}\delta b/(b_{0}-\delta b)$ [15]. We design the
functional form of $\delta b(y)$ such that the the capillary pressure
$p_{c}(x,y)$ within the defect is a linear function of $y$, i.e., $\delta
p_{c}(x,y)=\Pi(y-y_{d})$ for $y_{d}<y\leq y_{d}+\ell_{s}$, and $\delta
p_{c}(x,y)=\Pi\ell_{s}$ for $y_{d}+\ell_{s}<y<y_{d}+\ell$. The length of the
sloping part of the defect is $\ell_{s}=10$ mm and $\ell_{p}=\ell-\ell_{s}=5$
mm is the length of the constant capillary pressure zone (plateau) after the
slope. The defect width is $w=10$ mm. The value for the pressure slope
$\Pi=dp_{c}/dy$ was chosen through a series of experiments as $\Pi=6\times
10^{3}$ Pa/m. This value was chosen as it provides interface deformations that
are large enough to be easily captured by image analysis, while not too large
to avoid highly nonlinear deformations and snap-off events during drainage.
Under these conditions, the defect profile is: $\delta
b(y)=b_{0}-\left(\frac{\Pi(y-y_{d})}{2\gamma}+\frac{1}{b_{0}}\right)^{-1}$ for
$y_{d}<y\leq y_{d}+\ell_{s}$, and $\delta b(y)=\delta b(\ell_{s})$ for
$y_{d}+\ell_{s}<y<y_{d}+\ell_{s}+\ell_{p}$ (here we made the assumption that
the liquid perfectly wets the medium). The cell’s area is 6 cm $\times$ 6 cm,
with a depth of $b_{0}=3.6$ mm. For wetting fluid we used silicone oil, with
kinematic viscosity $\nu=10$ cSt $=10$ mm2/s, surface tension against air
$\gamma=20$ mN/m and density $\rho=0.93$ g/mL. The non-wetting fluid is
ambient air at atmospheric pressure.

The experiment is driven by changing the height of a reservoir of silicone oil
connected to the inlet (bottom) of the model. As the reservoir height
increases by $\delta h$, the oil pressure at the bottom of the model increases
by $\delta p=\rho g\sin(\alpha)\delta h$ thus causing the interface to move
upwards (imbibe), and vice versa in drainage. A moving average filter is
applied to the $h_{f}$ and $\eta_{m}$ data (in Figs. 7 and 9) to remove
spurious high-frequency noise (a consequence of image analysis inaccuracies).

The experiments in Fig. 7 _qualitatively_ demonstrate that (1) the imbibition-
drainage cycle around a single defect can be reversible, thus characterizing a
weak defect, and (2) the interaction between weak defects can trigger
irreversibility (hysteresis). The experiments also reveal further intriguing
physics not considered in our model. This is evident in the reversible case
shown in Fig. 9b: Examining the entire curve in Fig. 9b, outside the gray
shaded area which corresponds to the region in Fig. 7 where potential
hysteresis and dissipation may occur according to our model (see Fig. 2),
exposes other interesting features. The blue arrow in Fig. 9b points to a
perturbation in the curve even before the baseline position reaches the
defect. This can be explained through the 3-D nature of the liquid-air
interface. As the silicone oil wets both the top and bottom surfaces of the
cell, there is a thin film of liquid ahead of the 2-D projected interface.
Once this film touches the defect, the perturbation starts to grow even before
its 2-D projection touches the defect. Similarly, the red arrow in Fig. 9b
points to a small bump in the drainage cycle, which might be caused by pinning
of the contact line as it touches the defect in drainage. These 3-D effects
are not included in our 2-D model. The plateau region included after the
sloping part of the defect ensures that artifacts such as those related to
contact line pinning (red arrow in Fig. 9b) occur far from the area of
interest in our model (gray shading in Fig. 9b).

###### Acknowledgements.

RH acknowledges support from the Engineering and Physical Sciences Research
Council (EP/V050613/1); MD and JO received support from the Spanish Ministry
of Science and Innovation through the project HydroPore (PID2019-106887GB-C31
and C32, PID2022-137652NB-C41 and C42), and RP through the project
PID2021-122369NB-I00. RP and JO acknowledge AGAUR (Generalitat de Catalunya)
for financial support through project 2021-SGR-00450. MM acknowledges the
support from the Research Council of Norway through projects 262644 and
324555.

## References

  * Sahimi [2011] M. Sahimi,  _Flow and Transport in Porous Media and Fractured Rock_ (Wiley VCH Verlag GmbH, N.Y., 2011).
  * Albers [2014] B. Albers, Modeling the hysteretic behavior of the capillary pressure in partially saturated porous media: A review, Acta Mech. 225, 2163 (2014).
  * Haines [1930] W. B. Haines, Studies in the physical properties of soil. V. The hysteresis effect in capillary properties, and the modes of moisture distribution associated therewith, J. Agr. Sci. 20, 97 (1930).
  * Morrow [1970] N. Morrow, Physics and thermodynamics of capillary action in porous media, Ind. Eng. Chem. Res. 62, 32 (1970).
  * Berg  _et al._ [2013] S. Berg, H. Ott, S. A. Klapp, A. Schwing, R. Neiteler, N. Brussee, A. Makurat, L. Leu, F. Enzmann, J.-O. Schwarz, M. Kersten, S. Irvine, and M. Stampanoni, Real-time 3D imaging of Haines jumps in porous media flow, P. Natl. Acad. Sci. USA 110, 3755 (2013).
  * Cueto-Felgueroso and Juanes [2016] L. Cueto-Felgueroso and R. Juanes, A discrete-domain description of multiphase flow in porous media: Rugged energy landscapes and the origin of hysteresis, Geophys. Res. Lett. 43, 1615 (2016).
  * Helland  _et al._ [2021] J. O. Helland, E. Jettestuen, and H. A. Friis, A discrete-domain approach to three-phase hysteresis in porous media, Water Resour. Res. 57, e2021WR029560 (2021).
  * Måløy  _et al._ [2021] K. J. Måløy, M. Moura, A. Hansen, E. G. Flekkøy, and R. Toussaint, Burst dynamics, up-scaling and dissipation of slow drainage in porous media, Front. Phys. 9, 10.3389/fphy.2021.796019 (2021).
  * McClure  _et al._ [2021] J. E. McClure, S. Berg, and R. T. Armstrong, Thermodynamics of fluctuations based on time-and-space averages, Phys. Rev. E 104, 035106 (2021).
  * Primkulov  _et al._ [2020] B. K. Primkulov, J. Y. Y. Chui, A. A. Pahlavan, C. W. MacMinn, and R. Juanes, Characterizing dissipation in fluid-fluid displacement using constant-rate spontaneous imbibition, Phys. Rev. Lett. 125, 174503 (2020).
  * Bedeaux and Kjelstrup [2022] D. Bedeaux and S. Kjelstrup, Fluctuation-dissipation theorems for multiphase flow in porous media, Entropy 24 (2022).
  * Holtzman  _et al._ [2023] R. Holtzman, M. Dentz, R. Planet, and J. Ortín, The relation between dissipation and memory in two-fluid displacements in disordered media, Geophysical Research Letters 50, e2023GL104073 (2023).
  * Joanny and de Gennes [1984] J. F. Joanny and P.-G. de Gennes, A model for contact angle hysteresis, J. Chem. Phys. 81, 552 (1984).
  * Planet  _et al._ [2020] R. Planet, L. Díaz-Piola, and J. Ortín, Capillary jumps of fluid-fluid fronts across an elementary constriction in a model open fracture, Phys. Rev. Fluids 5, 044002 (2020).
  * Holtzman  _et al._ [2020] R. Holtzman, M. Dentz, R. Planet, and J. Ortín, The origin of hysteresis and memory of two-phase flow in disordered media, Commun. Phys. 3, 222 (2020).
  * Sethna  _et al._ [1993] J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhansl, B. W. Roberts, and J. D. Shore, Hysteresis and hierarchies: Dynamics of disorder-driven first-order phase transformations, Phys. Rev. Lett. 70, 3347 (1993).
  * Ortín and Goicoechea [1998] J. Ortín and J. Goicoechea, Dissipation in quasistatically driven disordered systems, Phys. Rev. B 58, 5628 (1998).
  * de Gennes [1986] P.-G. de Gennes, Imperfect Hele-Shaw cells, J. Physique (Paris) 47, 1541 (1986).
  * Soriano  _et al._ [2002] J. Soriano, J. Ramasco, M. Rodríguez, A. Hernández-Machado, and J. Ortín, Anomalous roughening of Hele-Shaw flows with quenched disorder, Phys. Rev. Lett. 89, 026102 (2002).
  * Geromichalos  _et al._ [2002] D. Geromichalos, F. Mugele, and S. Herminghaus, Nonlocal dynamics of spontaneous imbibition fronts, Phys. Rev. Lett. 89, 104503 (2002).
  * Darcy [1856] H. Darcy,  _Les fontaines publiques de la ville de Dijon._ (Victor Dalmont, Paris, 1856).
  * Saffman and Taylor [1958] P. G. Saffman and G. Taylor, The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid, Proc. R. Soc. A 245, 312 (1958).
  * Lavi  _et al._ [2023] I. Lavi, L. Rose, R. Planet, J. Casademunt, S. Santucci, and J. Ortín, Prebifurcation enhancement of imbibition-drainage hysteresis cycles, Phys. Rev. Fluids 8, 124002 (2023).
  * Sethna  _et al._ [2006] J. P. Sethna, K. A. Dahmen, and O. Perkovic, Chapter 2 – random-field ising models of hysteresis, in  _The Science of Hysteresis_, edited by G. Bertotti and I. D. Mayergoyz (Academic Press, 2006) pp. 107 – 179.
  * Middleton [1992] A. A. Middleton, Asymptotic uniqueness of the sliding state for charge-density waves, Physical Review Letters 68, 670 (1992).
  * Bertotti and Mayergoyz [2006] G. Bertotti and I. D. Mayergoyz, eds.,  _The Science of Hysteresis_ (Academic Press, 2006).
  * Goicoechea and Ortín [1994] J. Goicoechea and J. Ortín, Hysteresis and return-point memory in deterministic cellular automata, Phys. Rev. Lett. 72, 2203 (1994).
  * Keim  _et al._ [2019] N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry, and S. R. Nagel, Memory formation in matter, Rev. Mod. Phys. 91, 035002 (2019).
  * Bense and van Hecke [2021] H. Bense and M. van Hecke, Complex pathways and memory in compressed corrugated sheets, Proc. Nat. Acad. Sci. USA 118, 10.1073/pnas.2111436118 (2021).
  * Díaz-Piola  _et al._ [2017] L. Díaz-Piola, R. Planet, O. Campàs, J. Casademunt, and J. Ortín, Fluid front morphologies in gap-modulated hele-shaw cells, Phys. Rev. Fluids 2, 094006 (2017).
  * Note [1] We also note that our 2-D model is strictly valid only when the out-of-plane thickness varies slowly in space, while in the examples considered here the defects contain abrupt changes in thickness.
  * Lim and Saloma [2002] M. Lim and C. Saloma, Emergence of hysteresis in a network of nonhysteretic agents with continuous responses, Phys. Rev. Lett. 88, 038701 (2002).
  * Hunt  _et al._ [2013] A. G. Hunt, R. P. Ewing, and R. Horton, What’s wrong with soil physics?, Soil Sci. Soc. Am. J. 77, 1877 (2013).
  * Robbins and Joanny [1987] M. O. Robbins and J. F. Joanny, Contact angle hysteresis on random surfaces, EPL-Europhys. Lett. 3, 729 (1987).
  * Giacomello  _et al._ [2016] A. Giacomello, L. Schimmele, and S. Dietrich, Wetting hysteresis induced by nanodefects, Proc. Nat. Acad. Sci. USA 113, E262 (2016).
  * Lindeman and Nagel [2021] C. W. Lindeman and S. R. Nagel, Multiple memory formation in glassy landscapes, Sci. Adv. 7, eabg7133 (2021).
  * Wiese [2022] K. J. Wiese, Theory and experiments for disordered elastic manifolds, depinning, avalanches, and sandpiles, Rep. Prog. Phys. 85, 086502 (2022).
  * Shohat and Lahini [2023] D. Shohat and Y. Lahini, Dissipation indicates memory formation in driven disordered systems, Phys. Rev. Lett. 130, 048202 (2023).
  * Press  _et al._ [2007] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes_ , 3rd ed. (Cambridge University Press, Cambridge, UK, 2007).