# Stable and unstable capillary fingering in porous media with a gradient in grain size

Tom Vincent-Dospital, Marcel Moura, Renaud Toussaint, Knut Jørgen Måløy

► To cite this version:

Tom Vincent-Dospital, Marcel Moura, Renaud Toussaint, Knut Jørgen Måløy. Stable and unstable capillary fingering in porous media with a gradient in grain size. Communications Physics, 2022, 5 (1), pp.306. 10.1038/s42005-022-01072-1. hal-04311849

HAL Id: hal-04311849

https://hal.science/hal-04311849v1

Submitted on 28 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY-SA 4.0 - Attribution - ShareAlike - International
License
<figure><img src="image_2.png" /></figure>

# Stable and unstable capillary fingering in porous media with a gradient in grain size

Tom Vincent-Dospital<fnref n="1" />✉, Marcel Moura<fnref n="1" />✉, Renaud Toussaint<fnref n="1" />✉, <fnref n="1" />,²✉ & Knut Jørgen Måløy<fnref n="1" />✉, <fnref n="1" />,³✉

Multiphase flows in complex porous networks occur in many natural processes and engineering applications. We present an analytical, experimental and numerical investigation of slow drainage in porous media that exhibit a gradient in grain size. We show that the effect of such structural gradient is similar to that of an external force field on the obtained drainage patterns, when it either stabilises or destabilises the invasion front. For instance, gravity can enhance or reverse the drainage pattern in graded porous media. In particular, we show that the width of stable drainage fronts scales both with the spatial gradient of the necessary pressure for pore invasion and with the local distribution of this (disordered) threshold. The scaling exponent results from percolation theory and is – 0.57 for 2D systems. Overall, introducing a dimensionless Fluctuation number, we propose a unifying theory for the upscaling of dual immiscible fluid flows covering most classical scenarios.

<footnote n="1">¹PoreLab, The Njord Centre, Department of Physics, University of Oslo, Oslo, Norway. ²Université de Strasbourg, CNRS, Institut Terre & Environnement de Strasbourg, ENGEES, UMR 7063, Strasbourg, France. ³PoreLab, Department of Geoscience and Petroleum, Norwegian University of Science and Technology, Trondheim, Norway. ✉email: tom.vincent-dospital@fys.uio.no; marcel.moura@fys.uio.no; renaud.toussaint@unistra.fr; maloy@fys.uio.no</footnote>
Two-phase flow in porous media is important in a wide range of sectors, and holds applications in activities as diverse as food processing, management of underground water resources, sequestration of greenhouse gases or crude oil recovery. It is also crucial in everyday phenomena such as plant watering or making a cup of coffee. Owing to its broad impact but also to its complexity, it is a multidisciplinary subject that has been long investigated by physicists, geoscientists, hydrologists, chemists, biologists, and engineers. When one fluid displaces another immiscible fluid, the resulting structures in their distribution present various shapes and complexity$^{1–5}$, that can be compact, ramified and fractal$^{6,7}$. These structures are controlled by the forces, which actually drive the flow. These forces can for instance be viscous$^{1,4,5,8–10}$, capillary$^{1,9–12}$, gravitational$^{13–18}$ or a combination of these$^{19}$, and may depend on physical parameters such as the wetting properties of the solid-fluids system$^{3}$ or changes in the local geometry of the porous medium$^{20–23}$. In this work, we study the effect of the latter, that is, of the geometry of the solid matrix. In particular, we consider materials with an average gradient in their disordered porous structure. In nature, examples of such structured but disordered systems are graded geological bedding$^{24}$, for instance formed in rivers of varying stream intensity or by turbidite deposition, which are of interest to hydrologists and petroleum geoscientists. After a first conjecture that correlated structures in the solid matrix could help prevent flow instabilities$^{20}$, recent studies have considered the drainage of porous media presenting such gradients, first using regular matrices only displaying a structural gradient$^{21}$, and then matrices that additionally hold some quenched disorder$^{22}$. It has there been shown that the structural gradient could either stabilise or destabilise the drainage front. A generalised capillary number, taking into consideration the matrices’ structure, has then been introduced$^{21}$ to predict the effective stability of the flow. It has also been proposed that, if significant enough, the structural gradient could mitigate the capillary fingering arising from the quenched disorder$^{22}$.

We here follow a different theoretical approach based on percolation theory, and the generalisation of the recently reintroduced fluctuation number$^{10}$ to describe the effect of the noise on the width of the invasion front. We show that the gradient in grain size has a similar effect on the flow of the fluids than that of gravity (or, alternatively, any other external field). We illustrate this equivalence with drainage experiments in 3D printed transparent matrices. This technique allows both an accurate control of the experimental models’ structure and quenched disorder and the visualisation of the fluids’ distribution when a non-wetting phase displaces a wetting one in the pore space. Finally, we support the proposed theory with invasion percolation simulations reproducing the experimental results and the predicted scaling law between the width of the invasion front and the fluctuation number.

## Results

**Porous model and experimental set-up.** The experimental set-up and the 3D printed model for our drainage experiments are shown in Fig. 1. We describe, in the Methods section, the details of the set-up and the procedure we followed to print the model and run the flow experiments. This quasi two-dimensional porous matrix (of porosity ~70%) is composed of cylinders that are distributed in a monolayer using a Random Sequential Adsorption (RSA) algorithm$^{25}$. The RSA parameters are the minimum distance between the cylinders, the cylinders’ diameter, the cylinders’ height and the dimensions of the desired porous system. In our case, the RSA parameters vary between the model’s inlet and outlet. The total size of the model is 140 × 140 mm$^{2}$,

with the cylinder diameter $d$ varying linearly from 1 mm to 2 mm and the minimum cylinder separation varying from 0.4 mm to 0.8 mm, in order to preserve geometrical similarity along the model’s length. We thus design a model with a grain (cylinder) size gradient $\lambda \sim 0.007$ mm mm$^{-1}$. The height of the cylinder was chosen to be 2 mm. The model is then mounted in a flow cell, and, as per Fig. 1, the whole set-up can be tilted by an angle $\theta$ for the flow to occur in a chosen effective gravity field $g_0 \sin(\theta)$. Drainage experiments are then performed as described in the Methods section. A wetting water-glycerol solution, dyed with nigrosin, is displaced by air. The drainage is slow so that viscous forces are negligible, with a small capillary number $C_a = \mu V/\gamma$ that is about $10^{-8}$, where $V$ is the typical flow velocity, $\mu$ is the water-glycerol’s viscosity and $\gamma$ is the surface tension at the fluid-fluid interface. Note that, in the case where viscous forces would not be negligible, the gradient designed in the model structure could actually be seen, in a Darcy description, as a permeability gradient. In our case, however, it is rather a gradient in capillary threshold.

**Theory—gradient in pore throats and external fields.** The requirement for invasion into one pore neck by the non-wetting fluid is that the capillary pressure $p$ between the two fluids overcomes the capillary threshold value $\hat{p}_t$ of this pore neck

$$p(x, y) > \hat{p}_t(x, y), \qquad (1)$$

where $x$ and $y$ are the 2D coordinates for the pore neck position. In the case of our 3D printed model, the distribution in capillary pressure along the model is shown in Fig. 2. It is obtained from an approximation of Young-Laplace’s equation as $\hat{p}_t(x, y) = y \cos(\phi)(1/r + 1/h)$, where $r$ is the pore-throat radius, $h$ is half the cylinders’ height and $\phi$ is the contact angle – measured within the defending phase – for the solid matrix and the two fluids at play.

Assume that we have an external field that changes the capillary pressure linearly in the $x$ direction. One such field is the gravitational field (e.g., controlled by $\theta$ in Fig. 1). Another one, if the flow is fast enough, is the viscous pressure drop inside the fluid being withdrawn. In the case discussed in the present manuscript, the flow rate is slow enough for this viscous effect to be negligible compared to that of the other forces at stake. We will, however, keep a general formalism so they may be included. If we write the gradient in capillary pressure from the external fields as $G$, this capillary pressure at a position $(x, y)$ is

$$p(x, y) = p(x_0) + G \cdot (x - x_0), \qquad (2)$$

where $p(x_0)$ is the capillary pressure at an arbitrary position $x_0$. The condition for invasion can then be written as

$$p(x_0) > \hat{p}_t(x, y) - G \cdot (x - x_0) = p_t(x, y), \qquad (3)$$

where we have introduced the modified thresholds $p_t$. The system is now mapped onto a system without fields but where the thresholds are modified with the fields as a linear term in $x$. The capillary pressure $p = p(x_0)$ for this equivalent system is constant over the model. We now consider this modified system. The mapping between the occupation probability in percolation $f$ and the capillary pressure $p$ is given$^{26}$ by

$$f(x) - f_c = \int_{p_c(x)}^{p} N(p_t, x) dp_t \qquad (4)$$

In this expression, $f_c$ is the critical occupation probability, $p_c$ is the critical percolation pressure and $N$ is the normalised distribution in capillary pressure threshold. This distribution is a function of $x$ due to the overall gradient in pore neck and to the external field. Truly, such mapping between occupation probability and
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Experimental set-up. **a** Schematic cross-section of the experimental set-up. Two transparent PMMA plates (3 cm thick) are screwed together to close the 3D printed porous model with a gradient in grain and pore-throat sizes. Soft polymer layers (black layers) ensure good contacts between the PMMA plates and the model. The model is illuminated from below and pictures are taken from above. **b, c** Top-view pictures. Photographs of a non-wetting fluid (air) invading a wetting fluid (water-glycerol) in the square porous model of coordinate system (x,y). The side of the model is L = 140 mm. A pump is connected to the outlet and the inlet is open to air. The model can be tilted, leading to a capillary pressure gradient G = Δρg₀ sin(θ) in the x direction, where g₀ = 9.82 m s⁻² and Δρ ~ 1200 kg m⁻³ is the difference of density between the wetting and non-wetting fluids. The full-extent width of the invasion front highlighted in **c** and, to better characterise the average front width, we will denote by η the standard deviation of the front position along x.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Statistical description of the porous model. **a** Probability density function N̄ of ̂p_t(x, y)/̂p_c in the six areas (denoted 1 to 6) defined in Fig. 1 on the 3D printed model. This distribution is conserved in our model and, around ̂p_t = p_c it scales approximately as (̂p_t/̂p_c)⁻², as highlighted by the dash-dotted line. Note: the average pore connectivity in our model being 3, we defined the percolation pressure ̂p_c such that²⁶,³¹ ∫₀⁻¹ N(̂p_t) d̂p_t ≈ 0.65. **b** Evolution of the critical capillary pressure ̂p_c along the model (crosses), arising from its intrinsic structure. The straight line shows a linear fit, which is a first-order approximation of ̂p_c (the engineered gradient in the model was actually in grain size and not in percolation pressure). To this gradient, one can compute an equivalent model tilt (see Fig. 1): θₐₗ = arcsin(-∂p_c(x)/∂x / Δρg₀). Here, θₐₗ ~ 1.35 degrees.</figcaption></figure>

pressure arises from the ordinary percolation theory rather than from invasion percolation with cluster trapping, which is typically at stake in draining porous materials (e.g., see Fig. 1c where trapped fluid clusters are visible). Both models have actually been pointed out to belong to different universality classes²⁷, while, by contrast, invasion percolation without trapping shares the same class as ordinary percolation²⁷. In this work, we are interested in the invasion front and its stability. By definition, trapped wetting

clusters are not a part of this front, so that the different trapping rules do not affect its shape and statistical features but only the size of the trapped clusters that are left behind it. The standard tools of percolation theory should thus apply to describe the front. The match between the theory we now develop and the invasion percolation simulations with cluster trapping that we present in the next sections will confirm this hypothesis.

Let us now consider a stable invasion. In this case, a position x₁ exists were the capillary pressure is equal to the critical capillary pressure such that

$$p = p_c(x_1). \qquad (5)$$

We can also Taylor expand N(p_t, x) around the critical p_c(x₁) in Eq. (4) keeping only the lowest order in p - p_c(x₁) such that

$$f(x) - f_c = N(p_c(x_1), x)[p_c(x_1) - p_c(x)]. \qquad (6)$$

Expanding p_c(x) to first order in x - x₁

$$p_c(x) = p_c(x_1) + \left. \frac{\partial p_c(x)}{\partial x} \right|_{x_1} (x - x_1), \qquad (7)$$

we get

$$f(x) - f_c = -N(p_c(x_1), x) \left. \frac{\partial p_c(x)}{\partial x} \right|_{x_1} (x - x_1). \qquad (8)$$

Now, ∂p_c(x)/∂x will contain one term from the gradient in the critical capillary pressure ∂̂p_c(x)/∂x of the porous medium plus the term G, which is due to the external linear field. This gives

$$f(x) - f_c = aN(p_c(x_1), x) \left( G - \left. \frac{\partial \hat{p}_c(x)}{\partial x} \right|_{x_1} \right) \frac{x - x_1}{a}. \qquad (9)$$

Here, we have introduced a as the typical length of a pore. We now choose x such that η = |x₁ - x|, where η is the width of the front. We here write η = x₁ - x > 0, which corresponds to an invasion flow that progresses against the x direction if f(x) < f_c. The opposite convention for x could have, of course, also been chosen. Furthermore, we use Sapoval's assumption²⁸, that η scales in the same way as the correlation length ξ in percolation, (ξ/a) ∝ |f - f_c|⁻⁵, where v is a critical exponent typically equal to 4/3 for 2D systems²⁶. We obtain (η/a) ∝ (f_c - f)⁻⁵ and then, with Eq. (9), we find

$$\eta/a \propto F^{\frac{\nu}{1+\nu}}, \qquad (10)$$

where the exponent β = ν/(1 + ν) is about 0.57. For a 3D system,
$\beta$ would be approximately 0.47, with $\nu \sim 0.88^{26}$. We call the quantity $F$ the fluctuation number, which is given by

$$F(x) = a(x)N(p_c(x_1), x) \left( G - \left. \frac{\partial \hat{p}_c(x)}{\partial x} \right|_{x_1} \right). \quad (11)$$

$F$ is a dimensionless number dictating the invasion process. It is a generalisation of the fluctuation number introduced by Mâloy et al.¹⁰, Méheust et al.¹⁷ and Auradou et al.¹⁶. The quantity $1/N(p_c, x)$ characterises the typical width of the capillary threshold fluctuations, and, in the scenario of interest in our experiments, $a(G - \partial \hat{p}/\partial x)$ characterises the gravitational forces at the pore scale, corrected for the particular structure of the porous matrix. If the material disorder is important relatively to the gradient term, such a gradient becomes negligible at small scales (and reciprocally).

One can compare $F$ to other dimensionless numbers usually used to predict flow patterns in porous materials. In our case, a relevant one would for instance be the Bond number²⁹ $Bo \sim a^2 G/[\gamma \cos(\phi)]$, which compares the gravitational forces to the capillary ones. Previous works¹³,¹⁸ have proposed imbibition and drainage front widths to indeed scale as $Bo^{-v/(1+v)}$. Yet, and contrarily to $F$, the Bond number does not offer, a priori, any insight on the actual material disorder (e.g., $N$) or on an eventual structural trend (e.g., $\partial \hat{p}/\partial x$). We suggest that the Bond number could indeed sometimes help characterise the flow patterns in gravitational fields, because for regular enough geometries of the porous network (in particular without structural gradient), $Bo$ does scale as $F$ (see Supplementary Note 1). If one would, however, compare matrices with different distributions of disorder ($N$), such a direct scaling would fail. Furthermore, in the case where a matrix does hold a structural gradient, and where no external field applies (e.g., no gravity) the invasion can both be stable or unstable depending on the relative directions of the flow and the gradient, as we illustrate in the next sections. There, $Bo$ is - by contrast - always null, and an accurate stability prediction based on this number only would then be impossible.

As shown in Fig. 2a, our model was designed so that the dimensionless distribution $\tilde{N}(\hat{p}_t, \hat{p}_c(x)) = \hat{p}_c(x) \cdot N(\hat{p}_t, x)$ is conserved along the $x$ direction (i.e., it does not depend on $x$). Additionally, the pore throats in many media is to scale the same way as the size of the pores $a$. We can then write, as per the Young-Laplace law, $\hat{p}_c(x) \propto \gamma \cos(\phi)/a(x)$. Note that, in case of 2D models such as ours, this last expression is only strictly valid if the height of the system is big enough respective to the pore size (i.e., $h \gg a$), which in our model is correct for the smallest pores but is a brave assumption for the biggest ones for which $h \sim a$. This assumption, however, allows to rewrite Eqs. (10) and (11) as:

$$\eta \propto \left[ \frac{\tilde{N}(1)}{\gamma \cos(\phi)} \left( G - \left. \frac{\partial \hat{p}_c(x)}{\partial x} \right|_{x_1} \right) \right]^{-\beta} \times a^{-2\beta+1}. \quad (12)$$

Here, the dependence in $x_1$ of the gradient in percolation pressure has been removed, assuming a model where this pressure is perfectly linear in $x$ (that is, a model where the first-order expansion leading to Eq. (7) is exact). The exponent $-2\beta+1$ is close to 0 so that the effect of the pore size $a$ on the width of the front is small compared to that of the other terms.

This last expression is only a particular case of Eqs. (10) and (11), where the front width does not significantly evolve as the invasion progresses. In a more general porous material, the spatial distribution in pore size $a(x)$ and/or in pore invasion threshold $N(p_c, x)$ would matter.

**Experimental results.** From Eqs. (10) and (11), we see that a negative gradient in the critical capillary pressure threshold will

Flat models:

<figure><img src="image_5.png" /><figcaption>Tilted with angle $-2\theta_{eq}$ (gravity reversed stability):</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Drainage experiments in a 3d printed quasi two-dimensional model. In each photograph, air has invaded the model from the top of the pictures, replacing a water-glycerol solution coloured with nigrosin. All sides of the models are of dimension $L = 140$ mm. In **a**, $\partial \hat{p}_c(x)/\partial x < 0$, which corresponds to a decrease in 'permeability' as the flow progresses (the air invades from the big pores to the small ones). The invasion is stable. In **b**, $\partial \hat{p}_c(x)/\partial x > 0$, which corresponds to an increase in 'permeability' with the flow progression and an unstable processes. In both cases, the model is horizontal ($G = 0$). **c** and **d** show the same experiments in a chosen destabilising or stabilising gravity field (respectively) with tilt angle $-2\theta_{eq}$ (see Fig. 1 and 2). The stability of the flow is thus reversed.</figcaption></figure>

stabilise the front while a positive gradient will make it more unstable. In this sentence, the underlying convention is that the direction defining the gradient is opposing the general flow direction. This is illustrated in Fig. 3 that shows two drainage experiments in such configurations. In Fig. 3a we have a negative gradient in the critical capillary pressure threshold and the front is stabilised while in 3b the gradient is positive (the inlet and outlet having been swapped) and the front is unstable with a growing finger. In both cases the external (gravity) field is null. The stabilisation/destabilisation of the front is also reflected in the breakthrough time, which was about three times smaller in the case with the destabilising gradient (the externally imposed flow rate was the same for both experiments). Numerically analysing the images, we indeed measured that, in Fig. 3b, only 20% of the wetting phase is drained, while 65% is drained in Fig. 3a.

In Fig. 3c and d, we show the same experiments in a particular gravity field. The model is tilted with an angle $-2\theta_{eq} = -2 \arcsin(-\partial \hat{p}_c(x)/\partial x)$ as explained in Fig. 2, which effectively inverses the sign of the $(G - \partial \hat{p}_c(x)/\partial x)$ term in Eq. (11).

We should here restate that Eq. (11) only applies to the cases where the front is stable (i.e., to Fig. 3a and d) as front stability was an underlying hypothesis (to write Eq. (5)). As expected, for these two different experiments for which the fluctuation number $F$ is the same, the features of the invasion front are very similar. In the case of unstable fronts, we suggest that, the fluctuation
number could characterise the width of the invading fingers, rather than the width of the front. Indeed, in gravitational drainage experiments, such a finger width was indeed proposed¹⁵,³⁰ to scale with the Bond number, with the same exponent $\beta$, which we have here considered. In the case of the drainage of a matrix with a structural gradient (that is of interest here), this hypothesis will be further verified with the help of invasion percolation simulations.

## Invasion percolation simulations

**Stabilising gradient.** Fully validating Eq. (12) with experimental results remained a challenge. Indeed, because of the maximum size (~30 cm) of the models we could print, and because of the minimum distance between printed grains (a fraction of a millimetre) before they tended to unexpectedly merge, the range of gradient in capillary pressure that could be investigated was small. Owing to the theoretical scaling between such gradient and the front width (i.e., Eq. (12)), the range in obtainable $\eta$ was even smaller. Therefore, to verify our theoretical framework, we ran a series of invasion percolation simulations with cluster trapping, known to represent well capillary invasion processes.

These simulations are performed on square lattices (indexed along $x$ and $y$ directions), which abide by the same conditions that underlie Eq. (12). At a given $x$, the distribution $N$ in invasion threshold $\hat{p}_t$ for a pixel is random and uniformly distributed. Similarly to Fig. 2 the mean values of these distributions follow an arbitrary gradient along $x$ but the width of the distributions is such that $\tilde{N}(p_t/\hat{p}_c)$ is conserved along this same direction. Additionally, the size $a(x)$ of a pixel scales as $1/\hat{p}_c(x)$. Here, because our matrices' connectivity is 4, the percolation pressure $\hat{p}_c$ is such that$^{26,31} \int_0^{\hat{p}_c} N(\hat{p}_t) d\hat{p}_t = 0.5$. More details on the simulations can be found in the Methods section and the numerical code is available on a public repository³².

We ran ten simulations in the stable invasion domain and without external field ($\partial\hat{p}_c/\partial x < 0$ and $G=0$), varying the gradient in percolation pressure, but also the average size of the pixel $a$ in the model and the local width of the distribution in invasion threshold $1/\tilde{N}(1)$. We thus covered four decades of $\tilde{N}(1) \times \partial\hat{p}_c/\partial x < 0$ and two decades of $a$. Additionally, we ran six simulations where a linear external field ($G$) modifies the local fluid pressure, in a way that the invasion remains stable (i.e., $G - \partial\hat{p}_c/\partial x > 0$, although the respective sign of both terms may vary from one simulation to the other). For each simulation, we extracted the simulated front width $\eta$, once it reached a plateau after enough time steps. We defined this width as the standard deviation of the $x$ coordinate of the front between the invading phase and the main cluster of the defending phase (see Fig. 4).

Finally, with a least squares method, we fitted Eq. (12) to the obtained data in order to invert for the $\beta$ exponent. This procedure provided a good fit of the simulated data with a coefficient of determination $R^2 \sim 0.998$, and we found $\beta \sim 0.56 \pm 0.02$, where the accuracy of the fit is computed by letting $R^2$ vary by 5%. This value for $\beta$ is satisfyingly close to our theoretical prediction for percolation theory ($\beta = 0.57$).

**Destabilising gradient.** We also ran similar simulations, but in the unstable case, that is, with $G - \partial\hat{p}_c/\partial x < 0$. We varied the same parameters (i.e., $\partial\hat{p}_c/\partial x$, $\tilde{N}$, $a$ and $G$). In this unstable configuration, rather than characterising the width of a stable front, we characterised the average width $W$ of the growing invasion finger. We defined such a width as $W = A/l$, where $A$ is the area occupied by the finger (i.e., the area surrounded by the plain red lines in Fig. 5, in which the results are shown) and $l$ is the length of the finger along the $x$ direction. Comparing Figs. 4 and 5, one can

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Invasion percolation simulations—stabilising gradient. **a** to **c**: invasion percolation maps for the three simulations denoted in **d**, that have different stabilising gradient in invasion threshold. The invading non-wetting phase occupies the white pixels. The plain (red) line is the drainage front on which the width $\eta$ is computed. Length are in an arbitrary unit, denoted a.u. **d**: Width of the front as a function of the pore size $a$, of the gradient $\partial\hat{p}_c/\partial x$, of the external field $G$ and of the local distribution in capillary pressure $\tilde{N}$ (logarithmic scales). Each data point was computed on an independent simulation. Crosses represent simulations with only a structural gradient in pore threshold ($G=0$), whereas circles mark simulations also including an external field ($G$ varying from $-10^{-4}$ to $+10^0$ a.u. depending on the simulation). The grey plane is a fit of Eq. (12). The slopes of this plane along the two axes are given by $-\beta$ and by $1-2\beta$, respectively, and both indicate $\beta = 0.56 \pm 0.02$.</figcaption></figure>

notice how similar is the scaling of $W$ and $\eta$. Although Eq. (12) was only formally derived for stable fronts, we then analogously write for the width of the fingers in the unstable scenario:

$$W \propto \left[ \frac{\tilde{N}(1)}{y \cos(\phi)} \left( -G + \frac{\partial \hat{p}_c(x)}{\partial x} \right) \right]^{-\beta'} \times a^{-2\beta'+1}, \quad (13)$$

when the same assumptions for the matrix structure than those underlying Eq. (12) are respected. Fitting this expression to our simulation results, we obtained a good match ($R^2 \sim 0.992$) for $\beta' \sim 0.52 \pm 0.05$. This value is close to the value of $\beta$. As previously mentioned, a similar scaling of $W$ with respect to the Bond number (rather than to $F$) was reported in experimental observations of unstable drainage fingers growing in a gravitational field¹⁵,³⁰,³³. Our simulations suggest that, when possible, one should consider our defined fluctuation number instead.
<figure><img src="image_8.png" /></figure>

**Fig. 5** Invasion percolation simulations—destabilising gradient. **a** and **b:** Invasion percolation maps for the two simulations denoted in **c**, that have different destabilising gradient in invasion threshold. The invading non-wetting phase occupies the white pixels. The plain (red) line is the drainage front on which the width *W* is computed. Length are in arbitrary unit, denoted a.u.. **c** Width of the finger as a function of the pore size *a*, of the gradient ∂p̃<sub>c</sub>/∂x, of the external field *G*, and of the local distribution in capillary pressure N̄ (logarithmic scales). Each data point was computed on an independent simulation. Crosses represent simulations with only a structural gradient in pore threshold (*G* = 0), whereas circles mark simulations also including an external field, (*G* varying from −10<sup>-3</sup> to −10<sup>-1</sup> a.u. depending on the simulation). The grey plane is a fit of Eq. (13). The slopes of this plane along the two axes are given by −β'<sup>′</sup> and by 1 − 2β'<sup>′</sup>, respectively, and both indicate β'<sup>′</sup> = 0.52 ± 0.05.

Note that, either for the stable or the unstable drainage case,
the most general formalism that has been introduced in this
article, for the scaling law of the flow patterns, is given by Eq. (10)
rather than by Eq. (12) or Eq. (13). Figure 6 then displays the
dependence of either η/a or W/a with the fluctuation number F in
our numerical simulations, showing good agreement with the
theory.

**Conclusion**

In this paper, we have discussed the importance of capillary
fluctuations in porous media, as well as the characteristic length
scales in two-phase flow patterns set by the competition between
capillary fluctuations, external fields (e.g., gravitational or viscous
ones) and a gradient in the matrix percolation pressure.

<figure><img src="image_9.png" /><figcaption>**Fig. 6** Scaling of the flow patterns with the fluctuation number.</figcaption></figure>

Illustration of expression (10), which is here the most general formalism for the scaling of the size of the flow pattern with the fluctuation number, as defined by (11). The numerical results of Figs. 4 and 5 are shown (crosses and circles) and fitted with power laws (dotted line) for, respectively, the stable (a) and the unstable (b) drainage cases. Crosses represent simulations with only a structural gradient in pore threshold (G = 0), whereas circles mark simulations also including an external field.

In the case of fluid fronts that are stabilised by these fields and
porous media geometry, the fluctuation number *F*, which
describes the scaling of the front width *η*, was introduced. There,
the derived scaling exponents directly result from percolation
theory. When considering a viscous and gravitational field, the
theory describes well the scaling of the width of the fluid front
and the final saturation of the fluid left behind the invasion front
observed in laboratory experiments¹⁰. As shown here, it can also
predict the stabilisation and destabilisation of the front width in
such experiments depending on the sign of the spatial gradient in
the critical capillary pressure. Truly, more experiments are nee-
ded to conduct a quantitative experimental investigation of the
dependence of the scaling of the front width η when a structural
gradient is present, the challenge being to obtain a permeability
gradient varying over several decades in the laboratory. Standard
invasion percolation simulations have, however, here underlined
how reasonable the predicted scaling law is. In the case of a
destabilised flow, these simulations also enabled us to infer a
similar scaling law for the width of growing drainage fingers.

On a length scale smaller than η, the structure within the front
is generally fractal, while on a length scale larger than η, it is
homogeneous. The characteristic length scale η should thus be of
primary importance in defining a relevant representative ele-
mentary volume (REV) for an average Darcy description of the
two-phase flow problem³⁴. We suggest that characterising the
fluctuation number F, for instance from drilled core samples in
geological contexts, would help in this prediction of the front
width. Estimating the actual distribution in pore throats, needed
to derive N, can indeed be achieved with various experimental
techniques³⁵,³⁶. One could also characterise the relationship
between the fluids' pressures and saturations for systems with
structural gradients, by extending the approach of Moura et al.³¹
and Ayaz et al.³⁴. A good knowledge of this relationship is indeed
relevant for reservoir geophysicists and hydrologists.

In these geological applications, one should also consider the
usual repeating sequences of sediment layers with a gradient in
permeability in a given direction. This is for instance often
observed in fluvial or turbidite deposition²⁴. There, if the gradient
inside a single unit is stabilising, a change in layer would yet
correspond to a local but brutal destabilising effect (and reci-
procally). Predicting the flow behaviour over large distances
would thus likely require to take into consideration the typical
wave lengths of such layer repetitions, and compare them to the
typical length scale η of the invading pattern. Such a study would
be a natural continuation of the present work.
**Methods**

**Experimental method.** To generate our porous model, we have used Formlabs Form 3L printer³⁷, which employs a stereolithography 3D printing technology to produce objects in a transparent plastic material (Formlabs Clear Resin), composed of methacrylic acid esters and photoinitiators. This technique is based on the laser polymerisation of the resin. It allows to control the geometry of the porous network and in particular to fine-tune its grain and pore-throat distributions, for example by introducing a gradient in the pore sizes. In practice, cylinders are distributed in a monolayer using a Random Sequential Adsorption (RSA) algorithm²⁵. The cylinders are printed within a frame (which is part of the total print) closing the sides and bottom of the porous medium so that only the top of the cylinders are not sealed. Pipe connectors are also included (i.e., printed) at the inlet and outlet of this frame to facilitate the connection of the model in the experimental set-up. The spatial resolution of the printed models is about 0.1 mm. All the RSA parameters are specified in the manuscript.

The model is mounted in a flow cell, which was constructed around the print, in a way that optimises the visualisation of the pores. The 3D printed model is inserted between two layers of a soft polymer (polyvinyl chloride, PVC-2mm thick) whose main role is to efficiently seal the top of each cylinder. Around these soft layers two 3-cm thick PMMA plates are screwed together to confine the flow. On the side, the model is closed by 3D printed walls. The tubular inlet and the outlet of the model feed large channels along the whole model width so that the boundary conditions are the same along this direction. Each layer is transparent and the flow cell lies on a white light box to allow a good quality imaging with a reflex camera.

The invading fluid is air and the defending one is a mixture of 20% water and 80% glycerol, where the percentages relate to the total mass. In this mixture a nigrosin dye has been added (4 g per litre of water) to create an imaging contrast between the air and the liquid. The coloured mixture is the wetting phase, with a contact angle to the solid resin of 60 ± 10° in presence of air (measured optically on several droplets). The wettability of the soft PVC plate sealing the top of the printed cylinders was inferred to be similar. With a syringe pump, the wetting liquid mixture is withdrawn at a constant flow rate from one of the two model's ends while the air, connected to the atmospheric pressure of the laboratory, invades the model. We used a 0.3 ml/h flow rate.

**Numerical method.** The invasion percolation simulations are performed on square lattices (indexed along x and y directions), which abide by the same conditions that underlie Eq. (12). At the initial simulation stage, solely the first line of the matrix (the inlet) is invaded, so that the initial invasion front is a straight line. The pixels' connectivity is four, so that a pixel is connected to its top, bottom, right and left neighbours. On the side of the matrix, pixels have only three neighbours, and, in the corners, only two. At each time step, one new pixel is invaded. This pixel is the neighbour of the invading phase that has the lowest invasion threshold and which is still connected to the matrix outlet (i.e., it is not trapped). The invasion front is recomputed after the invasion of each pixel and the simulation stops when the front reaches the last line of the matrix.

**Data availability**

High-quality experimental photographs are available to the reader at the following location³², DOI: 10.5281/zenodo.7076135. Further requests should be addressed to Knut Jørgen Måløy (maløy@fys.uio.no).

**Code availability**

The numerical code for the invasion percolation simulations is available to the reader at the following location³², DOI: 10.5281/zenodo.7076135. Further requests should be addressed to Knut Jørgen Måløy (maløy@fys.uio.no).

Received: 28 April 2022; Accepted: 1 November 2022;

Published online: 29 November 2022

**References**

1. Lenormand, R., Touboul, E. & Zaarcone, C. Numerical models and experiments on immiscible displacement in porous media. *J. Fluid Mech.* **189**, 165 (1988).

2. Lenormand, R. Flow through porous media: limits of fractal pattern. *Proc. R. Soc. Lond.* **A 423**, 159 (1989).

3. Zhao, B., Mac Minn, C. W. & Juanes, R. Wettability control on multiphase flow in patterned microfluidics. *Proc. Natl Acad. Sci. USA* **113**, 10251 (2016).

4. Måløy, K. J., Feder, J. & Jøssang, T. Viscous fingering fractals in porous media. *Phys. Rev. Lett.* **55**, 2688- (1985).

5. Chen, J.-D. & Wilkinson, D. Pore-scale viscous fingering in porous media. *Phys. Rev. Lett.* **55**, 1892 (1985).

6. Mandelbrot, B. B. *The Fractal Geometry of Nature* (W. H. Freeman, 1982).

7. Feder, J. *Fractals* (Plenum,1988).

8. Weitz, D., Stokes, J., Ball, R. & Kushnick, A. Dynamic capillary pressure in porous media: Origin of the viscous-fingering length scale. *Phys. Rev. Lett.* **59**, 2967 (1987).

9. Lovoll, G., Méheyst, Y., Toussaint, R., Schmittbuhl, J. & Måløy, K. J. Growth activity during fingering in a porous hele shaw cell. *Phys. Rev. E* **70**, 026301 (2004).

10. Måløy, K. J., Moura, M., Hansen, A., Flekkøy, E. & Toussaint, R. Burst dynamics, upscaling and dissipation of slow drainage in porous media. *Front. Phys.* **9**, https://doi.org/10.3389/fphy.2021.796019 (2021).

11. Lenormand, R. & Zarcone, C. Invasion percolation in an etched network: Measurment of a fractal dimension. *Phys. Rev. Lett.* **54**, 2226 (1985).

12. Primkulov, B. K., Zhao, B., MacMinn, C. W. & Juanes, R. Avalanches in strong imbibition. *Commun. Phys.* **5**, 52 (2022).

13. Wilkinson, D. Percolation model of immiscible displacement in the presence of buoyancy forces. *Phys. Rev. A* **30**, 520 (1984).

14. Birovljev, A. et al. Gravity invasion percolation in two dimensions: experiment and simulation. *Phys. Rev. Lett.* **67**, 584 (1991).

15. Frette, V., Feder, J., Jøssang, T. & Meakin, P. Buoyancy-driven fluid migration in porous media. *Phys. Rev. Lett.* **68**, 3164 (1992).

16. Auradou, H., Måløy, K. J., Schmittbuhl, J., Hansen, A. & Bideau, D. Competition between correlated buoyancy and uncorrelated capillary effects during drainage. *Phys. Rev. E* **60**, 7224 (1999).

17. Méheyst, Y., Lovoll, G., Måløy, K. J. & Schmittbuhl, J. Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects. *Phys. Rev. E* **66**, 051603 (2002).

18. Breen, S. J., Pride, S. R., Masson, Y. & Manga, M. Stable drainage in a gravity field. *Adv. Water Resour.* **162**, 104150 (2022).

19. Toussaint, R. et al. Two-phase flow: Structure, upscaling, and consequences for macroscopic transport properties. *Vadose Zone J.* **11**, https://doi.org/10.2136/vzj2011.0123 (2012).

20. Yortsos, Y. C., Xu, B. & Salin, D. Delineation of microscale regimes of fully-developed drainage and implications for continuum models. *Comput. Geosci.* **5**, 257 (2001).

21. Rabbani, H. S. et al. Suppressing viscous fingering in structured porous media. *Proc. Natl Acad. Sci. USA* **115**, 4833 (2018).

22. Lu, N. B., Browne, C. A., Amchin, D. B., Nunes, J. K. & Datta, S. S. Controlling capillary fingering using pore size gradients in disordered media. *Phys. Rev. Fluids* **4**, 084303 (2019).

23. Niblett, D., Niasar, V. & Holmes, S. Enhancing the performance of fuel cell gas diffusion layers using ordered microstructural design. *J. Electrochem. Soc.* **167**, 013520 (2019).

24. Kuenen, P. H. Significant Features of Graded Bedding. *AAPG Bull.* **37**, 1044 (1953).

25. Hinrichsen, E. L., Feder, J. & Jøssang, T. Geometry of random sequential adsorption. *J. Stat. Phys.* **44**, 793 (1986).

26. Stauffer, D. & Aharony, A. *An Introduction to Percolation Theory* (Taylor & Francis, London, 1994).

27. Wilkinson, D. & Barsony, M. Monte carlo study of invasion percolation clusters in two and three dimensions. *J. Phys. A Math. Gen.* **17**, L129 (1984).

28. Sapoval, E., Rosso, M. & Gouyet, J. The fractal nature of a diffusion front and the relation to percolation. *J. Phys. Lett.* **46**, L149 (1985).

29. Hager, W. H. Wilfrid noel bond and the bond number. *J. Hydraulic Res.* **50**, 3 (2012).

30. Vasseur, G. et al. Flow regime associated with vertical secondary migration. *Mar. Pet. Geol.* **45**, 150 (2013).

31. Moura, M., Florentino, E., Måløy, K. J., Schafer, G. & Toussaint, R. Impact of sample geometry on the measurement of pressure-saturation curves: experiments and simulations. *Water Resour. Res.* **51**, https://doi.org/10.1002/2015WR017196 (2015).

32. Vincent-Dospital, T., Moura, M., Toussaint, R. & Måløy, K. J. Code and experimental images - Stable and unstable capillary fingering in porous media with a gradient in grain size. https://doi.org/10.5281/zenodo.7076135 (2022).

33. Wagner, G., Birovljev, A., Meakin, P., Feder, J. & Jøssang, T. Fragmentation and migration of invasion percolation clusters: experiments and simulations. *Phys. Rev. E* **55**, 7015 (1997).

34. Ayaz, M., Toussaint, R., Schafer, G. & Måløy, K. J. Gravitational and finite-size effects on pressure saturation curves during drainage. *Water Resour. Res.* **56**, https://doi.org/10.1029/2019WR026279 (2020).

35. Minagawa, H. et al. Characterization of sand sediment by pore size distribution and permeability using proton nuclear magnetic resonance measurement. *J. Geophys. Res. Solid Earth* **113**, https://doi.org/10.1029/2007JB005403 (2008).

36. León y León, C. A. New perspectives in mercury porosimetry. *Adv. Colloid Interface Sci.* **76-77**, 341 (1998).
37. Formlabs. *Technical Information, Form 3L*, Tech. Rep. (Formlabs, 2020). https://mediaformlabs.com/m/5f0ed4f707037528/original/-ENUS-Form-3L-Manual.pdf

**Acknowledgements**

We acknowledge the support of the University of Oslo, the Njord Center, and the Research Council of Norway through the PoreLab Center of Excellence (project number 262644) and the FlowConn Researcher Project for Young Talent (project number 324555). We also thank the University of Strasbourg and the IRP France-Norway D-FFRACT. We thank Per Arne Rikvold, Eirik Grude Flekkøy and Alex Hansen from PoreLab for insightful discussions. We also thank Mihailo Jankov for his support of the experimental work.

**Author contributions**

K.J.M. and T.V.D. proposed the theory developed in this article, T.V.D. and M.M. performed the experimental work, T.V.D. wrote and ran the numerical simulations, and R.T. advised in both the theory and the numerical implementation. T.V.D. and K.J.M. wrote the first version of the manuscript and all the authors agreed on the submitted version. Readers are welcome to comment and correspondence should be addressed to maloy@fys.uio.no.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s42005-022-01072-1.

**Correspondence** and requests for materials should be addressed to Tom Vincent-Dospital, Marcel Moura, Renaud Toussaint or Knut Jørgen Måløy.

**Peer review information** *Communications Physics* thanks Ruben Juanes and the other anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022