Effects of Pore-Scale Disorder on Fluid Displacement in Partially-Wettable Porous Media

We present a systematic, quantitative assessment of the impact of pore size disorder and its interplay with flow rates and wettability on immiscible displacement of a viscous fluid. Pore-scale simulations and micromodel experiments show that reducing disorder increases the displacement efficiency and compactness, minimizing the fluid-fluid interfacial area, through (i) trapping at low rates and (ii) viscous fingering at high rates. Increasing the wetting angle suppresses both trapping and fingering, hence reducing the sensitivity of the displacement to the underlying disorder. A modified capillary number Ca* that includes the impact of disorder λ on viscous forces (through pore connectivity) is direct related to λ, in par with previous works. Our findings bear important consequences on sweep efficiency and fluid mixing and reactions, which are key in applications such as microfluidics to carbon geosequestration, energy recovery, and soil aeration and remediation.


Introduction
In many of the challenges we face today as geoscientists, in particular in the context of water and energy resources, fluid invasion into a porous soil or sediment is a key process.Examples include hydrocarbon migration and recovery, methane venting from hydrate-bearing sediments, drying and wetting of soils, and carbon geosequestration 1 .Fluid displacement has also gained attention as a scientific paradigm for problems where competitive domain growth and nonlinear interface dynamics are influenced by spatial heterogeneity of the system properties, such as magnetic domains and biological films 2 .
The complex interplay between capillary, viscous, and gravitational forces, wettability effects, and the underlying heterogenous pore geometry, leads to ramified, preferential flow paths or "fingering".Predicting the emergent patterns is challenging, because of the sensitivity to pore-scale details and the large number of coupled mechanisms and governing parameters which vary over a wide range of values and scales 3,4 .Fingering leaves pockets of the original, defending fluid, increasing the fluid-fluid interfacial area, with important consequences on efficiency of processes relying on fluid withdrawal such as enhanced hydrocarbon recovery 5 , or on fluid mixing and reaction such as subsurface remediation 6 or CO 2 sequestration 7 .
The practical importance of fluid displacement motivated intensive research, that identified several important mechanisms governing the displacement efficiency.In disordered media, slow displacement of a wetting fluid (by a less wetting fluid, termed "drainage") leads to a patchy pattern with multiple islands of defending fluid, called capillary fingering (CF).The process of forming these islands is called trapping, as lack of fluid continuity between them inhibits their mobilization.Increasing the flow rate (hence the strength of viscous forces relative to capillary ones, characterized by the capillary number Ca = µ d v/γ) modifies the pattern depending on the fluids viscosity ratio, M = µ d /µ i , towards (i) viscous fingering (VF) for an unfavorable ratio, M > 1, or (ii) compact (CO) for a favorable ratio (M < 1) 8,9 .Here µ is the viscosity, v is the characteristic fluid velocity, and γ is the interfacial tension, where subscripts d and i denote defending and invading, respectively.
Another important property is wettability, namely the relative affinity of the fluids to the solid.Increasing affinity of the invading fluid (from drainage to imbibition, when the invading fluid is more wetting) smoothes the invasion front, even at high, unfavorable M 10 .At low rates, increasing wettability leads to a transition between CF and CO, whereas increasing Ca enhances viscous instabilities promoting VF, regardless of wettability 11,12 .We note that the frequent use of the term "(in)stability" to describe the pattern, and the interchange of "stability" with "efficiency" or "compactness".To avoid confusion between patterns (e.g.VF) and underlying mechanisms (e.g.hydrodynamic instabilities), we use "compact", "efficient" or "stable pattern".
Experimental and computational capabilities limit the number of mechanisms that can be modeled.Here, we focus on the following subset: capillarity, viscosity, wettability, and their coupling with random disorder in pore and/or particle sizes.
We exclude here gravity and buoyancy effects.We also focus on size disorder alone, excluding other types of heterogeneity including in wetting properties 13 or with spatial correlations such as layers or dual porosity media 14,15 .

Current understanding of disorder effects
Disorder in pore sizes is inherent in natural media such as soils and rocks 4 , and, to some extent, even in engineered devices due to limited fabrication resolution 16 .In some industrial applications, where homogenous fluid distribution and control over resulting patterns are key, low disorder is crucial and highly desirable 17 .We note that disorder as well as preferential pathways and heterogenous fluid distributions can be advantageous in applications where fluid mixing and reactions 6,7 are needed, and undesirable when efficient sweep of hydrocarbons 5 or contaminants 6 are crucial.As summarized below, while the importance of disorder is widely-recognized, a systematic exploration including quantitative analysis of samples of different, known disorder at variable wettability and rates is lacking.
The effect of disorder on the invasion morphology when either capillary or viscous forces are dominant seems to be univocally accepted.For slow drainage, high disorder was shown to promote trapping, leading to a transition from CO to CF [18][19][20][21][22][23][24] .For slow imbibition, quasi-static simulations show a transition from compact to faceted (FA) growth as disorder is decreased 20,25 .Ordered FA patterns were also observed experimentally, for imbibition in ordered samples (uniform pore sizes) at favorable M 9 .An interesting correlation between fluid invasion and magnetic systems was drawn regarding the interplay between disorder and the smoothing ("stabilizing") action of cooperative events 26 .At the other extremity of high Ca (VF regime), decreasing disorder produces snowflake-like patterns 22,27 with fewer, less tortuous fingers ("ordered dendritic", OD).We stress that FA and OD are only obtained for samples with nearly-uniform pore sizes, here resulting from uniformly sized particles ordered on a lattice.Such patterns would not appear in granular materials (including monosized bead packs) due to the large variability in apertures provided by packing disorder.
At intermediate Ca where both viscosity and capillarity are important, the qualitative paradigm ("conventional wisdom" 24 ) is that lowering disorder enhances the displacement efficiency 18,[21][22][23][24] .Few works derived a modified capillary number, Ca * , which includes the impact of disorder via the width of the capillary thresholds, neglecting its effect on the viscous forces.Quantitatively, this implies an inverse relationship between Ca * and disorder, namely that increasing disorder (at a given rate, Ca) produces a more CF-like pattern, increasing the transitional Ca between CF to VF, Ca CF/VF 18, 21-24 .We point to a subtle yet crucial inconsistency in the above: if increasing disorder increases Ca CF/VF , promoting a more CF-like pattern, sweep efficiency should increase (since VF is much less efficient that CF 18,[21][22][23][24]28 ), contrasting the observations stated above.
Here, we address this problemtic paradigm via a combination of pore-scale simulations, micromodel experiments and scaling analysis.To the best of our knowledge, our study is the first systematic, quantitative investigation of the impact of disorder and its coupling with rates and wettability.Our findings improve our understanding, generalizing the phase diagram of fluid invasion including rates, viscosity ratio, buoyancy and wettability effects 8,11,12,28 by adding an additional axis-disorder.

Problem setting and methodology
We consider radial displacement of a viscous fluid by an inviscid fluid (high, unfavorable ratio M) at a fixed rate Q in samples of different disorder in particle sizes.Details of the simulations and experiments are provided in Methods and Supplementary Material.Briefly, our 2-D medium is made of "particles" (cylindrical pillars) of diameters d drawn from a uniform distribution of width λ , d ∈ [1 − λ , 1 + λ ]d, providing well-defined pores and throats (of size r, with a triangular distribution), see Fig. 1a-c.Experimentally, we use a plastic micromodel initially saturated with glycerol, withdrawn from the cell's perimeter (outlet), allowing air to invade at its center.Time-lapse photography provides the pattern evolution (Fig. 1f).The large data set required for the quantitative analysis, including 240 runs on samples of different disorder, at various rates and wetting properties, was generated numerically, with the pore-scale model of Holtzman and Segre 12 .

Simulated patterns
Our simulations capture the plethora of invasion behaviors and emergent patterns as the disorder, rate and wettability are varied.For high disorder (λ =0.82), we reproduce the transition from CF to VF as the rate increases in drainage 8 (wetting angle θ < 90 • , θ measured through the defending fluid, see Fig. 1c), and the smoothing ("stabilizing") effect of increasing θ towards imbibition (θ > 90 • ), reducing trapping 11,12 (Fig. 2).We show that lowering the disorder λ enhances the displacement efficiency, by inhibiting trapping and formation of thin, tortuous fingers at low and high rates, respectively.While in highly-disordered systems the pattern is random (showing no dependence on the underlying structured lattice), reducing λ amplifies the impact of the underlying lattice, suppressing both trapping and fingering, eventually resulting in the following transitions: (i) from CF to CO in slow drainage 11 ; (ii) from CO to FA in slow imbibition 25 ; and (iii) from VF to OD 22,27 at high rates.We point that compact displacement is obtained despite of the highly unfavorable viscosity ratio, M = µ d /µ i ≈ 55.Reducing λ also suppresses the impact of wettability; that is, the competing effects of reducing λ and increasing θ offset each other (Fig. 2).

Experimental validation
To validate our numerical findings, we conducted micromodel experiments with two cells of different disorder (λ =0.22 and 0.52) of geometry identical to the simulations, displacing glycerol by air ( The experimental patterns support our findings that reducing λ enhances the displacement efficiency by (i) suppressing trapping at low Ca, leading to a transition from CF to CO; and (ii) reducing the number of tortuous fingers, increasing their width (from VF to OD) at high Ca (Fig. 3).The experiments also show the impact of λ on the invasion dynamics: (i) the transition in the mode of front propagation from radially-symmetric and continuous (CO) to intermittent and asymmetrical (CF) as λ is increased at low rates (Videos 1 and 2 in the the supplementary material); and (ii) radial growth of viscous fingers with vanishing impact of the underlying lattice as λ increases (from OD to VF, see Videos 3 and 4) at high rates.
We note that while patterns in simulations and experiments are similar in type (e.g.number and width of fingers or interface roughness), they are not identical.Potential causes for this mismatch include (1) manufacturing defects in pillar placement and shape (Fig. 1f, insets) and (2) mechanisms not included in our model such as snapoff 9,31 or droplet fragmentation 29 in low-porosity rocks, and thin wetting fluid films progressing ahead of the main imbibition front 30 .Manufacturing errors could greatly affect the pattern due to its sensitivity to small geometrical details; e.g.enlarging a single aperture can allow invasion of multiple pores, otherwise inaccessible.In contrast, for the current settings of relatively high viscosity contrast M, low contact angle (vs.150 • in 30 ), high porosity and small throat to pore size ratio, the aforementioned mechanisms should not be dominant 9,31 .A more conclusive validation requires improved manufacturing accuracy as well as a more rigorous, quantitative analysis.Such experiments are underway.Modeling mechanisms such as droplet fragmentation or film flow is a challenge which is outside the scope of this paper 3,4,31 , however of great interest for future research.

Pattern characteristics
We characterize the patterns through their finger width w, interfacial areas A inter and A front , and breakthrough saturation S. The interfacial areas (lengths in 2-D) are computed from the ratio of interfacial and invaded pores, including trapped regions (A inter ) or excluding them (A front ).The areas are normalized by the invaded volume V inv (area in 2-D), and multiplied by the lattice spacing a to make them dimensionless.The finger width w is evaluated by skeletonizing the pattern with a Voronoi algorithm, and measuring the distance (in lattice units a) of the medial line from the interface.In evaluating w we include small clusters trapped within fingers, unlike some previous studies where these were ignored 28 (resulting in larger w).
Our simulations point to the persistent effect (across variable Ca and θ ) of increasing disorder λ , manifested as an increase in A inter and A front (approaching the theoretical limit of 1 for VF at high rates), and a decrease in w (Fig. 4a-c).The counteracting effects of increasing θ and lowering λ are demonstrated by a nearly-constant, low A inter and A front and high w for λ < 0.6, in particular at low Ca (approaching the limit a/L for CO), see Fig. 4d-f.As shown below, this insensitivity is caused by two competing effects of decreasing λ in imbibition (high θ ): (i) reducing the variation in capillary thresholds, increasing compactness; and (ii) inhibiting cooperative events which smooth the interface, hence decreasing compactness.

Invasion selectivity and cooperativity
To further explain the coupling between disorder, rate and wettability, we define two characteristics: invasion selectivity and cooperativity.Selectivity, evaluated here from the mean invaded throat size r inv , increases with λ and Ca, while decreasing with 4/11 θ to a minimum at slow imbibition (Fig. 5a).Reduced selectivity implies greater dependency of the invasion pattern on the underlying pore geometry (e.g.FA and OD in Fig. 2).This occurs in ordered systems (λ =0.22), where the radially-symmetric invasion (Video 1 in supplementary material) samples nearly-equally all throat sizes (retaining the triangular distribution, Fig. 5c) such that the mean invaded size is close to the mean (r=0.54a).In contrast, high disorder provides constrictions with high entry pressures that remain non-invaded (Video 2 in supplementary material), leading to a skewed distribution with a larger mean r inv .Such high selectivity is mostly pronounced at slow drainage, where small pores are completely avoided (Fig. 5a-b).In imbibition, since cooperative invasion occurs at lower capillary pressures than burst instabilities 12 , smaller pores can be invaded than in drainage, reducing r inv .
Cooperativity ε coop is the relative number of nonlocal, cooperative pore filling events (overlaps).Overlaps, which increase in occurrence with θ (and do not appear below ∼60 • ), stabilize imbibition patterns by suppressing fingering and trapping 12 .Schematic side panels in Fig. 6 show two pairs of adjacent menisci at a similar curvature for λ =0.82 and 0.22, demonstrate that pore-size homogeneity hinders overlaps.Since disorder also promotes trapping, the two effects counteract, leading to the small sensitivity to λ in slow imbibition; this is particularly noticeable for A front , which is more affected by trapping than A inter or w (Fig. 4d-f).Increasing rates diminishes both selectivity and cooperativity, as pressure screening results in thin fingers invading only into more permeable regions 12 (Figs.5-6).Consequently, the mechanism for invasion selectivity changes with rates: at low Ca it is the high capillary thresholds along the interface (locally), whereas at high Ca the flow resistance-a nonlocal effect which depends on the connectivity of larger pores-becomes dominant.

Trapping and sweep efficiency
One of the most important characteristics of fluid displacement from a practical standpoint-describing the ability to produce oil or water and displace a contaminant-is the sweep efficiency.Sweep is intimately related to trapping, leaving behind isolated patches of the defending fluid.We define sweep S and trapping fraction χ trap from the following dimensionless volumetric ratios: S = V inv /V tot , where V tot is the volume of the convex polygon wrapping the pattern (Fig. 7a, inset), and χ trap = V trap /V inv , where V trap is the volume of trapped fluid (Fig. 7c-f, in red).
The complex interplay between disorder, wettability and rates is demonstrated by the nonlinear, non-monotonic behavior of sweep and trapping 32 .Both lowering disorder λ and increasing wetting angle θ provide a similar effect: suppressing trapping hence increasing sweep, to almost 100% at low rates (Fig. 7a-b).Flow rate strongly impacts both S and χ trap .In slow drainage (low Ca and θ ) and high λ , the patchy CF-like pattern traps a significant portion of the defending fluid (Fig. 7b; see pattern in   In imbibition, a significant portion of the pores (measured by the cooperativity ε coop ) are invaded by cooperative pore filling events ("overlaps") 12 .Lowering disorder hinders cooperativity, as demonstrated by the side panels showing two adjacent menisci with similar curvatures.At high flow rates, increasing dominance of viscous instabilities suppresses overlaps 12 , decreasing ε coop .Lines and shading show ensemble mean and standard deviation among four realizations.Fig. 7c).At slow imbibition (high θ ), lowering disorder inhibits both cooperative pore filling (Fig. 6) and trapping, effects which offset each other, explaining the observed minimal impact of λ (cf.Fig. 4d-f).Increasing the rate promotes viscous instabilities and invasion by thin fingers which leave most of the defending fluid behind (Fig. 7f), minimizing sweep despite of the reduced trapping.Rates change not only the trapped volume, but also the underlying mechanism: from multiple patches within thick fingers (CF) at low Ca (Fig. 7c-d) to fewer patches trapped between thin fingers (VF) at high Ca (Fig. 7f).

Scaling analysis: Introducing disorder effects
In light of the complex effect of disorder, and, in particular, the ambiguity in the literature regarding its interplay with viscous forces, we address this analytically via scaling, deriving a disorder-dependent number Ca * .We begin with the classical definition, Ca = δ p vis /δ p cap , where δ p vis and δ p cap are the characteristic viscous and capillary forces 9 .The force driving the growth of individual fingers perpendicular to the interface (along the direction of the externally-applied pressure drop) δ p vis is evaluated from the pressure drop in the viscous defending fluid over a characteristic length L vis , δ p vis ∼ ∇p vis L vis .Here ∇p vis ∼ µ d v/k is the gradient and k is the medium's permeability.The capillary force is evaluated from a characteristic entry pressure, δ p cap ∼ γ/r c , where r c is a characteristic aperture, providing 33  The classical expression 9 Ca = µ d v/γ is obtained by evaluating δ p vis at the pore-scale, L vis ∼ r, together with the simplifying assumptions r c ∼ r and k ∼ r 2 .

Ca
Here, we point to an inherent shortcoming in the above definition when disorder is introduced: the assumption k ∼ r 2 ignores the effect of interpore connectivity, which controls the resistance to viscous pressure dissipation (flow).We claim that disorder in pore throat sizes λ affects both (i) capillary forces δ p cap through the entry threshold distribution which controls the invasion order along the interface (locally), as well as (ii) viscous forces through the permeability k, a nonlocal effect, not accounted for in previous works 18,[21][22][23][24] .To include the latter in our definition of Ca * , we evaluate k as an effective permeability: considering finger growth in the radial direction, we use the weighted harmonic mean of the local (interpore) permeabilities k i ∼ r i 2 (i.e.resistors in series), k −1 ∼ 1/r 2 , providing The main outcome of Eq. ( 2) is that since increasing disorder (the width of r distribution) decreases the size of the smallest throats ("bottlenecks") and thus the permeability k (which is biased towards the smallest r values), it increases δ p vis .If the accompanying increase in δ p cap is smaller, this would amount to increasing Ca * .In our derivation, unlike the classical one presented above, the terms r c and r in Eq. ( 2) represent different characteristic lengths (though both at the pore-scale), hence are not equal: r c arises from the characteristic threshold, δ p cap ∼ γ/r c , whereas r is our choice for the pore scale characteristic length for pressure dissipation, L vis ∼ r.Furthermore, the choice of r c in Eq. ( 2) depends on wettability.In drainage, since invasion occurs mostly though burst instabilities 12 which strongly depend on throat sizes r, the distribution is dominant; thus, we evaluate r c using the standard deviation, r c = σ (r).In imbibition, the increasing occurrence of overlaps with λ (Fig. 6) makes for a complex dependence of the characteristic invaded size r c on disorder.Since predicting this dependence is out of the scope of this paper, we simply use r c = r for imbibition.
As samples with higher λ exhibit more tortuous, thinner fingers (at a given Ca, see insets in Fig. 8), plotting the pattern characteristic A inter vs. our rescaled Ca * partially collapses curves of different λ , (Fig. 8; similar collapse observed for A front and w).This indicates a direct relationship between Ca * and λ (Fig. 8b, lower inset), namely a lower transitional Ca CF/VF value at higher λ , in par with the inverse relationship (and increase in Ca CF/VF ) in the literature 18,[21][22][23][24] .Our scaling analysis also quantifies the decreased sensitivity of the emergent pattern to λ as the wetting angle θ increases (highlighted in Fig. 4d-f): Ca * /Ca is ∼10 times larger for λ =0.82 than for λ =0.22 for drainage (θ =5 • ), vs.only ∼2.5 for imbibition (θ =120 • ), see inset in Fig. 8b.We note that while our scaling is a valuable quantitative explanation of our findings, the expression in Eq. ( 2) and consequent values (Fig. 8) were derived for our 2-D medium with uniform distribution of round particles; extension to other types of media is of interest for future research.We explain quantitatively the impact of disorder through scaling analysis.Plotting interfacial area A inter against our rescaled capillary number Ca * (rather than Ca, insets) provides a partial collapse of curves of different λ .This indicates a direct relationship between Ca * and λ (b, lower inset), suggesting a stronger impact of disorder on the viscous forces, thus a lower transitional Ca CF/VF at higher λ .Our scaling also demonstrates quantitatively the diminishing effect of λ as the wetting angle θ increases (b, lower inset; see text).Lines and shading depict the ensemble mean and standard deviation.

Discussion
We study the impact of random pore size disorder on the fluid-fluid displacement in porous media.We use a novel pore-scale model which provides a mechanistic description of both drainage and imbibition.Our numerical results are validated against micromodel experiments, and explained theoretically by scaling analysis.To the best of our knowledge this is the first systematic, quantitative study of the impact of disorder and its interplay with flow rates and wettability.
We show that lowering disorder λ increases the pattern compactness, with a smoother fluid-fluid interface.We demonstrate that the effect of disorder depends on its coupling with flow rate and wettability: (i) at low rates, a more pronounced impact of λ on capillary thresholds and hence on the order of invasion and trapping; (ii) at high rates, increasing λ strongly increases the viscous pressure drop by decreasing pore connectivity, promoting pressure screening and viscous instabilities, increasing the number of thin, tortuous fingers; and (iii) increasing the wetting angle θ stabilizes the invasion pattern by reducing trapping and smoothing the interface, thus minimizing the impact of disorder.At low λ , the underlying pore geometry (lattice structure) becomes dominant, in particular in imbibition, leading to symmetrical growth of faceted, compact pattern at low Ca and dendritic fingers at high Ca.
While our observation that decreasing disorder increases the sweep efficiency agrees with previous studies 18,[20][21][22][23][24][25]27 , our interpretation does not: we emphasize the previously ignored impact of pore-scale disorder on the effective, sample-scale connectivity (permeability), and thus on the viscous forces. Wequantify this effect by deriving a disorder-dependent capillary number Ca * providing a direct relationship between Ca * and λ , that is a lower transitional Ca CF/VF value at higher λ , in par with previous works 18,[21][22][23][24] .Our findings delineate the complex interplay among disorder, rates and wettability, via a simplified 2-D model.Valuable future work includes generalization of the model towards 3-D granular pack with correlated heterogeneities (more representative of soils and rocks), inclusion of additional mechanisms such as droplet fragmentation or film flow, and a more rigorous experimental validation by improving manufacturing accuracy and quantitative analysis.
In summary, our work offers a fresh perspective into a long-standing paradigm of significance in processes where efficiency of fluid extraction/injection and mixing and reaction of the fluids are crucial.Increasing disorder promotes fingering, which limits the recovery of subsurface fluids such as hydrocarbons 5 or contaminants 6 , while increasing interfacial areas which controls fluid mixing and reactions 34 in microfluidics 16,17 , subsurface remediation 6 , soil aeration 35 , and long-term CO 2 geosequestration by carbonate mineralization 7 .

Model medium: 2-D analog
The main objective of this work is to gain fundamental understanding of the complex interactions among disorder, wettability and rates.Therefore, we choose the simplest representation (for both numerics and experiments) of a partially-wettable, disordered porous medium that retains the most crucial aspects of the behavior of such medium.We consider here the following model system: cylindrical solid pillars placed on a regular lattice (Fig. 1), noting that consideration of the convoluted 3-D pore geometry (relevant for instance to natural sediments) would hinder basic insights, and thus should be left to a later stage.We introduce disorder through variations in particle sizes, rather than through randomness in particle positions or in the coordination number of pores, since the former allows a more rigorous quantitative comparison among different samples.Similarly, our choice of triangular lattice limits the number of menisci which can invade into a pore simultaneously to two, retaining the crucial aspect of cooperative pore filling while avoiding the ambiguity of resolving a burst on a square lattice 20 .
We construct square-shaped samples of size L = N x a by placing N x × N y particles (N y × 2N x pores), where N y =2N x / √ 3 and a is the lattice spacing.Each particle triplet defines a pore with volume V connected to three neighboring pores by throats with aperture 0 < r a; to avoid closed throats due to particle overlap we enforce d < a/(1 + λ ), see Fig. 1c.Pillar diameters are much smaller than their height, d h, such that the in-plane curvature is much larger than the out-of-plane curvature, R in R out .This allows us to approximate the meniscus as 2-D, with a radius of curvature R

Numerical simulations
Simulations are conducted with the model of Holtzman and Segre 12 , summarized below (see supplementary material for further details).This model can accommodate a wide range of properties and conditions including flow rates, wetting properties and viscosities.To the best of our knowledge, it is the only pore-scale model which provides a mechanistic description of the dynamics of partially-wetting invasion in a large, disordered domain.Previous models were limited by computational costs 3,4,31 , compromising the dynamics by assuming quasi-static displacement 19,20 or the size and heterogeneity of the domain 24,32 .We represent the fluid-fluid interface as a sequence of menisci, where each meniscus is a circular arcs which intersects a pair of particles at the prescribed contact angle θ , and has a radius of curvature R = γ/∆p according to the Young-Laplace law.Here ∆p is the pressure jump across the meniscus.The angle θ is measured through the defending fluid (i.e.θ < 90 • for drainage), and represents an effective advancing angle, which includes dynamic effects 3 .Knowledge of R and θ allows to resolve analytically the position and hence stability of each meniscus.
A crucial feature of our model is its account for cooperative pore-filling (overlaps 20 ), which become dominant at high θ .Overlap is a nonlocal mechanism as it is affected by menisci in multiple pores, acting to smooth the interface 12,20 .Menisci are tested for three types of capillary instabilities 20 : (1) Haines jump or burst, when the curvature exceeds a threshold set by the local geometry; (2) touch, when a meniscus intersects a third, downstream particle; and (3) overlap of adjacent menisci, destabilizing each other [Fig.1(c-e)].Once destabilized, the meniscus incipiently invades the downstream pore, at a rate evaluated from the fluid's viscous resistance (see supplementary material for details).The finite pore filling time in our model captures the invasion dynamics, including the spatiotemporal nonlocality associated with rapid interfacial jumps and the consequent interface readjustments and pressure screening 33,36 , overcoming a long-standing computational challenge 3,31 .

Experimental setup
Micromodels allow visualization of the solids and fluids, making them a highly useful tool to study fluid displacement.Here, we use 3-D printing ("stereolithography"), an additive method fabricating free-form solids using a computer-aided design (CAD) file.A ProJet 3500 Multi Jet Printer (3D Systems Inc.) with "VisiJet M3 Crystal" plastic provided in-plane and out-of-plane resolution of 50 and 32 µm, respectively 37 .The limited resolution resulted in rough, non-round pillars, and, though rare, a few that are misplaced (off-lattice), see insets in Fig. 1f.We fabricated two cells of geometry similar to the numerical samples (L = 100a) with λ = 0.22 and 0.52, however with lattice spacing a = 1 mm (vs.0.5 mm in simulations).Pillar height of h = 3 mm was chosen to ensure sufficiently large aspect ratio h/r such that the meniscus can be approximated as 2-D, while small enough to minimize hydrostatic pressure and thus 3-D, out-of-plane effects inducing vertical flows.The lattice is embedded in a frame with a trench (of width 1 cm) around it to ensure evenly distributed withdrawal.The chip is sandwiched between two 1 cm-thick acrylic plates, with an outlet cross-sectional area A out ≈ 5.5•10 −4 m 2 .
The cell is initially filled with glycerol (∼99.5% purity with µ d ≈ 1Pa•s and γ ≈ 6.3 • 10 −2 N/m at ∼25 • C, Gadot Chemical Terminals LTD.).A high-viscosity fluid was selected to probe high Ca, where the CF-VF transition occurs.The contact angle (for both acrylic and VisiJet M3 Crystal) was measured to be 75 • ± 5 • using a sessile drop goniometer (EasyDrop DSA20E, Krüss GmbH).Glycerol is withdrawn from ports at the cell corners by a syringe pump (Harvard Apparatus PHD ULTRA 70-3007) at a fixed rate Q, letting air (µ i = 1.8 • 10 −5 Pa•s) invade freely through a central inlet (Fig. 1f).Ca values in Fig. 3 correspond to Q=0.2, 0.5, 2, 5, 20, 50 and 100 ml/min.The cell is held horizontally to avoid gravity effects.The pattern evolution is tracked by time-lapse photography using a CMOS camera (uEye 3520, IDS 1.3MP).We illuminate the cell from above with LEDs.To improve image contrast, the bottom acrylic plate was painted black, and the inner side of the top plate was sandpapered, causing dry areas to appear opaque.

Figure 1 .
Figure 1.Our porous medium is made of variably-sized cylindrical pillars placed on a triangular lattice.(a) Numerically, we track the fluid-fluid interface (black line) and fluid pressures (increasing from blue to red).(b) Zoom in showing the lattice (spacing a) of particles (diameter d) and pores interconnected by throats (of width r, see panel c).The fluid-fluid interface is represented as a sequence of circular menisci, touching particles at contact angle θ , with radius of curvature R set by the local capillary pressure.Menisci can be destabilized by: (c) burst; (d) touch; or (e) overlap.Brown arrows indicate direction of advancement, destabilized arc in dash (From Holtzman and Segre 12 ; Copyright 2015 by the American Physical Society).(f) Experimental setup: A micromodel is initially saturated with glycerol, which is withdrawn by a syringe pump from the cell's perimeter (outlet), allowing air to invade at its center (inlet).Time-lapse images provide the pattern evolution.Limited fabrication resolution results in pillars with rough, non-round edges, and, in a few places, misplaced (inset, in red; a = 1 mm).

Figure 2 .
Figure2.Simulated displacement patterns at various flow rates (Ca) and wetting properties (θ ), for three degrees of disorder, λ .Increasing λ decreases the displacement efficiency, in a manner that depends on both Ca and θ : promoting trapping and capillary fingering (CF) at slow drainage (low Ca and θ ), and viscous fingering (VF) at high rates (regardless of θ ).Reducing λ enhances the effect of the underlying geometry (the ordered lattice), which, together with the stabilizing effect of wettability results in faceted (FA) and ordered dendritic (OD) patterns at low and high Ca, respectively.Pores filled with invading and defending fluid are marked in black and white, respectively (solid not shown).Sample size is L=260a.

Figure 3 .
Figure 3. Experimental displacement patterns patterns in micromodels with different disorder λ , showing that decreasing λ provides a more ordered pattern, with a transition from CF to CO and from VF to OD at low and high rates, respectively, in agreement with our numerical simulations.In these difference images, showing changes since air invasion started, air appears black, and solid pillars and glycerol (the defending fluid) are invisible (white).

Figure 4 .
Figure 4. Simulated pattern characteristics: (a) fluid-fluid interfacial area A inter ; (b) front area (excluding trapped regions) A front ; and (c) mean finger width w.The improved displacement efficiency as λ is lowered is manifested by a decrease in A inter and A front , and increase in w.Efficiency is also enhanced by decreasing rate (Ca) and increasing wettability (θ ).Slow imbibition (Ca=3.6•10−5 , θ =120 • , light blue in d-f) produces a compact front which is nearly unaffected by disorder.Lines and shading in a-c depict the ensemble mean and standard deviation among four realizations (with L=100a).Panels d-f show data from eight samples (λ ) and three rates (darker and thicker lines indicate higher Ca).

Figure 5 .
Figure 5. Increasing disorder increases the invasion selectivity, namely the mean invaded pore throat size r inv (panel a) by avoiding small pores (b).At low λ , the invasion samples nearly equally all throat sizes (c), such that r inv ≈ r = 0.54a.Selectivity is further reduced with the wetting angle θ , due to increasing occurrence of overlaps which are less sensitive to aperture sizes than bursts.In contrast, selectivity increases with rates, as screening produces thin fingers which completely bypass low permeability regions.Lines and shading in (a) are ensemble mean and standard deviation among four realizations.Panels b and c show the probability density function (pdf) of invaded throat sizes (r inv , colored lines), and the triangular throat size distribution (r, black dots) for two samples, λ =0.82(b) and 0.22 (c).

Figure 6 .
Figure 6.In imbibition, a significant portion of the pores (measured by the cooperativity ε coop ) are invaded by cooperative pore filling events ("overlaps")12 .Lowering disorder hinders cooperativity, as demonstrated by the side panels showing two adjacent menisci with similar curvatures.At high flow rates, increasing dominance of viscous instabilities suppresses overlaps12 , decreasing ε coop .Lines and shading show ensemble mean and standard deviation among four realizations.

Figure 7 .
Figure 7.Both sweep efficiency S and trapping fraction χ trap depend on the flow rate Ca, wetting angle θ and disorder λ in a complex, non-monotonic manner.Sweep efficiency S improves with decreasing λ and Ca and increasing θ (a).Rate has the largest impact on sweep, as viscous instabilities leave most of the defending fluid in place, greatly reducing sweep despite the reduced trapping, χ trap (b).The trapping mechanism changes with Ca from multiple patches distributed within the thick, capillary fingers (c-d) to fewer patches in between the thin viscous fingers (f); both mechanisms are suppressed by increasing θ and decreasing λ .Lines and shading in a-b depict the ensemble mean and standard deviation.Patterns in panels c-f (trapped and invaded pores in red and gray) refer to points on the curve for λ = 0.82 and θ =5 • in b (blue solid line).

Figure 8 .
Figure 8.We explain quantitatively the impact of disorder through scaling analysis.Plotting interfacial area A inter against our rescaled capillary number Ca * (rather than Ca, insets) provides a partial collapse of curves of different λ .This indicates a direct relationship between Ca * and λ (b, lower inset), suggesting a stronger impact of disorder on the viscous forces, thus a lower transitional Ca CF/VF at higher λ .Our scaling also demonstrates quantitatively the diminishing effect of λ as the wetting angle θ increases (b, lower inset; see text).Lines and shading depict the ensemble mean and standard deviation.