# Burst Dynamics, Upscaling and Dissipation of Slow Drainage in Porous Media

Knut Jørgen Måløy, Marcel Moura, Alex Hansen, Eirik Grude Flekkøy, Renaud Toussaint

► To cite this version:

Knut Jørgen Måløy, Marcel Moura, Alex Hansen, Eirik Grude Flekkøy, Renaud Toussaint. Burst Dynamics, Upscaling and Dissipation of Slow Drainage in Porous Media. Frontiers in Physics, 2021, 9, pp.796019. 10.3389/fphy.2021.796019. hal-03547298v2

HAL Id: hal-03547298

https://hal.science/hal-03547298v2

Submitted on 28 Jan 2022

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# Burst Dynamics, Upscaling and Dissipation of Slow Drainage in Porous Media

Knut Jørgen Måløy¹,²*, Marcel Moura¹, Alex Hansen³, Eirik Grude Flekkøy¹,⁴ and Renaud Toussaint⁵,¹

¹PoreLab, Department of Physics, The Njord Centre, University of Oslo, Oslo, Norway, ²PoreLab, Department of Geoscience and Petroleum, Norwegian University of Science and Technology, Trondheim, Norway, ³PoreLab, Department of Physics, Norwegian University of Science and Technology, Trondheim, Norway, ⁴PoreLab, Department of Chemistry, Norwegian University of Science and Technology, Trondheim, Norway, ⁵CNRS, Institut Terre et Environnement de Strasbourg, UMR 7063, University of Strasbourg, Strasbourg, France

We present a theoretical and experimental investigation of drainage in porous media. The study is limited to stabilized fluid fronts at moderate injection rates, but it takes into account capillary, viscous, and gravitational forces. In the theoretical framework presented, the work applied on the system, the energy dissipation, the final saturation and the width of the stabilized fluid front can all be calculated if we know the dimensionless fluctuation number, the wetting properties, the surface tension between the fluids, the fractal dimensions of the invading structure and its boundary, and the exponent describing the divergence of the correlation length in percolation. Furthermore, our theoretical description explains how the Haines jumps' local activity and dissipation relate to dissipation on larger scales.

**Keywords:** porous media, burst dynamics, intermittency, drainage, multiphase flow, critical phenomena, percolation

## 1 INTRODUCTION

Two-phase flow in porous media is crucial in a variety of sectors, ranging from fundamental research to applications in a wide array of industrial sectors such as fuel cell [1] and solar cell technology [2], fiber-reinforced composite materials [3], textile fabric characterization [4], prospection and exploration of oil and gas [5–8] etc. The description of flows inside natural porous media, such as soils and rocks, is also crucial for the study of groundwater flows [9, 10] and the treatment of soil contaminants [11–13] but it also matters for everyday tasks such as making a cup of coffee [14]. It is a multidisciplinary subject that has been investigated for decades by hydrologists, physicists, chemists, geoscientists, biologists, and engineers due to its practical importance and complexity. The structures observed are controlled by the forces involved, such as viscous [15–21], capillary [18–25], and gravitational forces [26–31], as well as wetting properties [32–38] and changes in the local geometry of the porous medium [39, 40]. The structures vary in shape and complexity [15, 16, 18, 32, 37, 41–47], from compact to ramified and fractal [48, 49]. The fractal nature of porous media is itself important for a number of applications [50, 51], such as electrolyte diffusion through charged media [52] a topic of relevance for the development of modern battery technology [53].

In most practical applications of porous medium physics, the typical length scales where our interest lies is substantially larger than the scale where the relevant physics is taking place. Oil reservoirs or water aquifers are in the range of kilometers while the typical pore sizes are commonly in the micrometer scale, about nine orders of magnitude smaller. How do we deduce the flow behavior at large scales from small-scale physics? The usual way of solving these problems

**OPEN ACCESS**

**Edited by:**

Fernando A. Oliveira,
University of Brasília, Brazil

**Reviewed by:**

Boqi Xiao,
Wuhan Institute of Technology, China
Hans Herrmann,
UMR7636 Laboratoire de physique et
mécanique des milieux hétérogènes
(PMMH), France

*Correspondence:
Knut Jørgen Måløy
k.j.maløy@fys.uio.no

**Specialty section:**
This article was submitted to
Interdisciplinary Physics,
a section of the journal
Frontiers in Physics

**Received:** 15 October 2021
**Accepted:** 15 November 2021
**Published:** 10 December 2021

**Citation:**
Måløy KJ, Moura M, Hansen A,
Flekkøy EG and Toussaint R (2021)
Burst Dynamics, Upscaling and
Dissipation of Slow Drainage in
Porous Media.
Front. Phys. 9:796019
doi: 10.3389/fphy.2021.796019
is a top-down approach using Darcy's law [54] on a mesoscopic level. However this approach does not take into account local fluctuations, like capillary or viscous fluctuations, which are averaged out. In this manuscript we take an alternative bottom-up approach where we emphasize on the pore-level capillary fluctuations and compare those fluctuations with the characteristic forces which are set up by the external fields on the whole system. Examples of such forces are gravitational or viscous fields. We will also limit our discussion to the stable drainage regime, which occurs when a nonwetting fluid displaces a wetting fluid, and when the viscous and/or gravitational forces stabilize the displacement front between the two fluids. This approach was introduced in the late 1970s [55, 56] and its theoretical development benefited greatly from invasion percolation models [57]. On this subject, other relevant experimental, numerical, and theoretical articles have since been published [26, 27, 58–61]. The dimensionless fluctuation number F, introduced in [30, 61], quantifies the ratio between the viscous and/or gravitational field and the capillary pressure fluctuations. The characteristic length scale η which describes the width of the invasion front, and which depends on F, is of central importance to calculate the saturation behind the front, which in turn gives a measure of the sweeping efficiency of a given drainage process, a quantity of great interest in a number of applications. The local structures will be fractal on length scales smaller than η with a crossover to a homogeneous behaviour on length scales larger than η. Knowing the scaling of the fluid structures up to the length scale η allows one to have full control of the energy balance of the problem and calculate large-scale quantities such as the final saturation behind the front, the dissipation or the entropy production, and the work required to move the front forward.

The dynamics of slow displacement in porous media has been observed to occur in an intermittent manner by so called Haines jumps [59, 62–72]. The invasion percolation model [55–57] describes well the structure of slow displacement in porous media [22], but it does not describe the dynamics realistically because the invasion is limited to one pore at a time. However it is possible to introduce a realistic interpretation of time in a modified invasion percolation model [63] by introducing a constant $\kappa$ relating the volume change at the interface with a corresponding change in pressure. References [63, 64] found both in simulations and experiments that the pressure fluctuations of the Haines jumps follow a power law distribution with an exponential cutoff function.

In the present work we use a new approach to calculate the dissipation and the energy of the surfaces left behind the front by considering the elastic energy (surface energy) released by the Haines Jumps when the front is in a steady state regime. In between jumps all the external work goes into building up the elastic (surface) energy of the fluid front. When the invading fluid front is in a steady state, its average length is constant, and the work applied to the system (for example by an external pump) must be equal to the elastic energy released by the invasion front (green line in Figure 1). As a result, using the distribution function of the capillary pressure fluctuations and performing a bottom-up approach to integrate up the elastic energy released by the bursts, we can check the consistency of our theory as well as

the distribution of dissipative events. We found as expected that the total elastic energy released by the bursts is equal to the work $W = \langle p \rangle \Delta V$, where $\langle p \rangle$ is the average pressure across the model and $\Delta V$ is the volume change of the invading fluid corresponding to the interface motion. This result provides an important consistency check for the analysis. The dissipation can then be calculated by subtracting the generated surface energy from the total applied work. The surface energy is directly computed by the scaling of the invasion structure, which is given by the invasion structure's fractal dimension, also measured in our work. As a result of the theory and quasi-two-dimensional experiments, we discovered that the work $W$, dissipation $\Phi$, and the energy of the surfaces left behind the front $E_s$ are all proportional to the number of pores $S$ invaded. Up to the characteristic length scale $\eta$, the spatial scaling of $S$ is fractal and for larger length scales the scaling becomes proportional to the system size instead. We then show

<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Upper figure: A nonwetting fluid 1 invades another wetting fluid 2 in a porous model of width $w$, length $L$ and coordinate system $(x,y)$. A syringe pump is connected to the lower side and the model is open to air on the upper side. The model can be tilted with an effective gravitational constant $g = g_0 \sin(\theta)$ along the x direction of the model, where $g_0 = 9.82 \text{ m/s}^2$. The average position of the front is $h$ and the width of the front is $2\eta$. Lower figure: Cross section of the experimental model system (A) shows two transparent Plexiglass plates, (B) the 3D printed porous model, (C) a PVC film kept under pressure to close the model and (D) screws that clamp the model. The model is illuminated from below and pictures are taken from above.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | The figure illustrates a nonwetting fluid invading a wetting fluid at a single pore-throat at different radius of curvature $r_0$ and $r = r_0 - dr$. The distance moved between the two front positions is $dx_i$ and $\psi$ is the wetting angle.</figcaption></figure>

that we can explicitly calculate $W$, $E_s$, and $\Phi$ if we know the surface tension of the fluids and the wetting properties from the experiments. We further discovered an important analytical result: that the ratios $E_s/W$ and $\Phi/W$ are both independent of system size.

## 2 EXPERIMENTAL TECHNIQUE

The majority of the experimental results presented in this paper are taken from data published in the past in our group. The system is composed of a modified Hele-Shaw cell [73], filled with a monolayer of glass beads with diameters in the range (1.0 mm, 1.2 mm). The glass beads are randomly distributed in the cell gap and the voids in between the beads form the porous network. This quasi-two-dimensional geometry allows for the direct visualization of the fluid phases, by means of regular optical imaging using a digital camera. For details about the model construction, see for instance Refs. [19, 74].

In this paper we also present results produced using an alternative stereolithography 3D printing technique to create the porous structures. We have employed a Formlabs Form 3 printer to produce models in a transparent plastic material (Clear Resin FLGPCL04). This technique allows us to control the geometry of the porous network and in particular to tune its porosity. In the experiments presented here we have made quasi-two-dimensional models where cylinders are distributed with a Random Sequential Adsorption algorithm [75] where we can set the minimum distance between the cylinders (see Figure 1), typically chosen as 0.3 mm. The cylinders height and diameter were both chosen to be 1 mm. The spatial resolution of these models is about 0.09 mm and the models are constructed to optimize the visualization of the pores, which are seen from a top-down view. The 3D printed porous model is placed between two thick Plexiglass plates which are clamped around the edges using screws to give robustness to the setup and ensure the quasi-two-dimensional geometry. A flexible PVC film is placed between the porous model and the top Plexiglass plate. This film plays an important role: due to its flexibility, when the screws around the

<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | The figure illustrates the dependence of the pressure across the model $p(t)$ and the time $t$. The front is in equilibrium between point 1 and point 2 and reaches the threshold value of one of the pores along the front at point 2. A burst will appear between point 2 and 3 and the front is there out of equilibrium.</figcaption></figure>

model are fastened the PVC film gets in contact with the top of all cylinders, thus ensuring the appropriate sealing of the model. The film has similar wetting properties as the 3D printing material used in the construction of the model. Figure 1 shows a typical snapshot of an experiment and a diagram of the setup. We also define in the upper part of the figure the model's width $w$, length $L$, the invading front (green line), its average position $h$ and width $2\eta$.

The porous network is initially fully saturated with a wetting viscous liquid composed of a mixture of glycerol (80% in weight) and water (20% in weight). The kinematic viscosity, density and surface tension (with respect to air) are $\nu = 4.25 \cdot 10^{-5} \text{ m}^2/\text{s}$, $\rho = 1.205 \text{ g/cm}^3$ and $\gamma = 0.064 \text{ N/m}$. The wetting liquid is dyed with a dark blue colorant (Ligroin), to aid visualization. Air is used as the nonwetting phase. The contact angle measured inside the wetting phase is $\psi = 70^\circ$. During an experiment, the liquid phase is withdrawn with a syringe pump (Harvard Apparatus) at a constant flow rate, leaving the model from a width-spanning channel at the bottom end of the cell. Air enters the model from the top, through another width-spanning channel that is open to the atmosphere.

## 3 SURFACE ENERGY, DISSIPATION AND BURST DYNAMICS

Consider a single pore identified by the index $i$ in which a nonwetting fluid displaces a wetting fluid as illustrated in Figure 2. The front moves from a position with a radius of curvature $r_0$ to one with a radius of curvature $r = r_0 - dr$. This results in a volume change $dv_i$ of the nonwetting fluid, and a corresponding volume change—$dv_i$ of the wetting fluid in that pore. Assume that we measure the distances in the normal (radial) direction from each of the points on the front defined with radius $r_0$ to the front defined by the radius of curvature $r =$
$r_0—dr$. The longest of these distances we denote as $dx_i$. The distance $dx_i$ to lowest order in $dr$ is given by

$$dx_i = \left. \frac{dx_i}{dr} \right|_{r_0} dr = \alpha_i dr, \qquad (1)$$

where

$$\alpha_i (\psi, p_0) = \left. \frac{dx_i}{dr} \right|_{r_0} \qquad (2)$$

depends on the wetting angle $\psi$ (See **Figure 2**) and the capillary pressure $p_0$ given by the Young-Laplace equation $p = \gamma(1/r + 1/r_1)$ [10] for $r = r_0$. Here $r$ is the in-plane, and $r_1$ the out of plane radius of curvature assumed to be constant in the quasi-two-dimensional experiments.

We can then calculate the corresponding increase in capillary pressure

$$dp = \frac{\gamma dr}{r_0^2}. \qquad (3)$$

Then by using Eqs 1, 3 we get

$$dp = \frac{\gamma}{r_0^2 \alpha_i} dx_i. \qquad (4)$$

The volume change $dv_i$ of the invading fluid in the pore can be written as

$$dv_i = A_i (\psi, p_0) dx_i, \qquad (5)$$

where $A_i (\psi, p_0)$ is the surface area in pore $i$ which is function of the wetting angle and the capillary pressure $p_0$. We then get

$$dp = \frac{\gamma}{r_0^2 A_i \alpha_i} dv_i, \qquad (6)$$

such that

$$dv_i = \kappa_i dp, \qquad (7)$$

where the capacitive volume

$$\kappa_i (\psi, p_0) = \frac{r_0^2 A_i (\psi, p_0) \alpha_i (\psi, p_0)}{\gamma}, \qquad (8)$$

gives the fluid volumetric change per unit capillary pressure in pore $i$ [63].

Let $U_i(x)$ be the elastic energy (surface energy) of the interface in pore $i$, and $dx_i$ a small displacement of the interface between the two fluids due to an external pump driving the system at a low flow rate such that viscous forces can be neglected. Assume that the system is in mechanical equilibrium. The work performed by the pump will then increase the elastic energy due to this displacement. The elastic energy associated to pore $i$ having interface located at the position $x_i$ would then be

$$U_i (x_i) = U_i (0) + \int_0^{x_i} p A_i dx_i, \qquad (9)$$

where $p$ is the capillary pressure across the interface, and $U_i(0)$ is the elastic energy when $x_i = 0$. Using Eqs 7, 9, and 5 we get

$$U_i (x_i) = U_i (0) + \kappa_i \int_{p_0}^{p} p^* d p^* = U_i (0) + \frac{\kappa_i}{2} p^2. \qquad (10)$$

We have assumed $p_0 = 0$ without loss of generality since $\kappa_i p_0^2/2$ also can be included in $U_i(0)$. If instead of considering a single pore we consider the work done on a front having $n$ pores, the total elastic energy $U_t$ will be the sum of the contributions from all pores belonging to the invasion front (green line in **Figure 1**). We then get

$$U_t = U_t (0) + n \frac{\kappa}{2} p^2, \qquad (11)$$

where

$$\kappa = \left( \sum_{i=1}^{n} \kappa_i \right) / n, \qquad (12)$$

where $\kappa_i$ is the capacitive volume of pore $i$ so that $\kappa$ is the average capacitive volume over the $n$ interface pores. The capillary pressure $p_2$ at the time $t_2$ and the capillary pressure $p_1$ at the time $t_1$ are related by

$$p_2 - p_1 = \frac{q(t_2 - t_1)}{n\kappa}, \qquad (13)$$

where $q$ is the imposed volumetric flow rate. Since we have assumed the system to be in mechanical equilibrium, the interface will slowly increase its capillary pressure without any bursts (Haines jumps) [62]. The capillary pressure will then build up from $p_1$ at $t_1$ to $p_2$ at $t_2$ due to the work $W$ performed by the external pump, see **Figure 3**. This work will then increase the total elastic energy (surface energy) of the interface from $U_t$ ($t_1$) to $U_t$ ($t_2$)

$$U_t (t_2) - U_t (t_1) = W = n\kappa \int_{p_1}^{p_2} p^* d p^* = \frac{n\kappa}{2} (p_2^2 - p_1^2). \qquad (14)$$

After some time the interface will however reach a situation where the capillary pressure in one of the pore-throats is at the threshold value $p_t$. The interface at that pore will then become unstable, and the invading fluid will move into one or more neighboring pores. At this time, the capillary pressure at the interface of the growing burst will be lower than the capillary pressure at the other parts of the interface. This produces a local velocity field that extends from the growing burst to the other parts of the interface. The interface will then back-contract, beginning with the pores closest to the pore where the burst begins and spreading out through the interface until it reaches equilibrium. Then the capillary pressure is again the same at all pores along the interface. Let us assume that the burst starts at time $t_2$ at a capillary pressure $p_2 = p_t$ and that, after the burst, the system reaches another equilibrium state at capillary pressure $p_3$. During the burst the elastic energy of the part of the fluid interface which is back-contracting will be reduced. This energy reduction is equal to
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**FIGURE 4** | (A) Tracking of the invasion front (green) and of the in-plane invasion contours left behind the invasion front (red). Notice that some isolated beads are included in the tracking. (B) Measurement of mass fractal dimension *D*, the fractal dimension *D*<sub>e</sub> of the invasion contours and the fractal dimension *D*<sub>f</sub> of the main invasion front (green line). The dotted and dashed lines have slopes of -2 and -1 respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**FIGURE 5** | Tracking of the internal interfacial contours between the nonwetting and wetting phases (orange) and the nonwetting and solid phases (green).</figcaption></figure>

$$U_t(p_2) - U_t(p_3) = \frac{\kappa n(t_2)}{2} (p_2^2 - p_3^2) = \\ \frac{\kappa n(t_2)(p_2 - p_3)(p_2 + p_3)}{2}, \qquad (15)$$

and will go to creation of the surface energy of the new burst in addition to viscous dissipation due to the flow that takes place during the burst. Here $\Delta p = p_2 - p_3$ is the pressure drop in the burst and $(p_2 + p_3)/2$ is the average pressure during the burst.

The surface energy of one single burst in our quasi-two-dimensional experiments has two contributions, one from the area of the burst $2\gamma_2 a^2 s$ and one from the interface contour $a\gamma_1 b s_e$. Here *a* is the average in-plane pore size, *b* is the height of the

cylinders (or the diameter of the beads in the case of the previously published experiments based on a monolayer of glass beads, see **Section 2**), *s* is the number of pores in the burst considered, and $a s_e$ is the length of the invasion contour of the burst (including trapped clusters). The factor 2 in the term $2\gamma_2 a^2 s$ is due to the creation of two new solid-air interfaces after the burst, one at the bottom and one at the top. The need for two different surface energies values, $\gamma_1$ and $\gamma_2$, is explained below.

In **Figure 4A**) we show the tracking of the invasion front (green line), and the in-plane invasion contour left behind the front (red line). The fractal dimensions *D* of the invading fluid phase (air), $D_e$ of the invasion contours, and $D_f$ of the invasion front are shown in **Figure 4B**), where a box counting technique was employed in the measurement [49]. We have obtained the values $D = 1.87 \pm 0.10$, $D_e = 1.84 \pm 0.10$ and $D_f = 1.45 \pm 0.10$. The dotted and dashed lines on the left correspond to the exponents -2 and -1 respectively. These are expected as the fractal nature of the invasion structure has a lower bound at the pore-size *a*. For boxes of size $\delta < a$, we recover the intrinsic 2 dimensional nature for the invading structure mass (dotted line) and the 1 dimensional nature of the perimeters of the front and internal contours (dashed lines).

When a fluid film is left behind the invasion front the surface tension $\gamma_1 = \gamma_2 = \gamma$ which is the surface tension between the two fluids. However when no film is left behind the invasion front, $\gamma_2 = \gamma_{ns} - \gamma_{ws}$, where $\gamma_{ns}$ is the interface tension between the nonwetting fluid and the solid and $\gamma_{ws}$ is the interface tension between the wetting fluid and the solid. Using the Young equation (10) we can also write the previous relation as $\gamma_2 = \gamma \cos(\psi)$, where $\psi$ is the static contact angle at the liquid-solid-air triple line. In the situation in which no wetting fluid film is left behind the invasion, the in-plane interface contour (red lines in **Figure 4**)
<figure><img src="image_9.png" /><figcaption>**FIGURE 6** | Burst size distribution G(s) from slow drainage experiments. The line has a slope of -1.32 corresponding to the results predicted by the analytical prediction Eq. 21. The data for EXP-2 (blue triangles) has been shifted vertically by 0.4. Data from Ref. [71].</figcaption></figure>

will be partly formed by nonwetting-wetting segments and partly formed by nonwetting-solid phase segments. This division is exemplified in Figure 5, where we have split the contours shown in red in Figure 4 into orange and green segments, where orange denotes the interface between the nonwetting and wetting phases and green the interface between nonwetting and solid phases. Let us define the ratio $\epsilon = A_{ns}/A_t$, where $A_{ns}$ is the interface's surface area between the nonwetting phase and the solid phase, and $A_t$ is the total surface area of the interface. In Figure 5, $\epsilon$ corresponds to the ratio between the total length of green lines to the total length of both green and orange lines added. We then get an effective surface tension $\gamma_1 = (1 - \epsilon)\gamma + \epsilon(\gamma_{ns} - \gamma_{ws})$ for the total internal interfacial contours. Using again the Young equation (10) we can rewrite this expression as $\gamma_1 = (1 - \epsilon)\gamma + \epsilon\gamma\cos(\psi)$.

We will in the following consider a situation where we have a long channel containing the porous medium with the invading fluid entering from one of the short sides of this channel and with the outlet on the other side. We will further assume that we follow the displacement front. The front will have a width of the order of the width $w$ of the channel. Apart from an initial transient we can then consider the front to be in a steady state regime with a constant average length $n$ (measured in number of pores). Since the elastic energy of the growing interface on average must be constant, the work $W$ applied by the pump must equal the elastic energy released by the bursts which is equal to the dissipation $\phi$ plus the surface energy needed to create the new interface $E_s$, such that

$$W = \Phi + E_s. \qquad (16)$$

It is important to note that $W$, $\Phi$, and $E_s$ are average quantities that do not account for fluctuations, and that $\phi$ and $E_s$ are averaged over a large number of bursts.

From [63] we know that the distribution function of the burst sizes $G(s)$ is directly linked to the distribution of the capillary pressure drops $\mathcal{F}(\Delta p)$ during the bursts since

$$s = n\kappa\Delta p/\bar{v}, \qquad (17)$$

<figure><img src="image_10.png" /><figcaption>**FIGURE 7** | The dependence of the crossover function $H(\Delta p/\Delta p^*)$ of the pressure jump distribution in a modified invasion percolation simulation [63]. For each value of $\kappa$, the points represent averages over five independent simulations on 200 × 1,500 lattices. The scaling exponent 0.30 = 1/(1 + $\nu D$). Data is taken from [63].</figcaption></figure>

where $\bar{v}$ is the average single pore volume averaged over the invading structure, and

$$\mathcal{F}(\Delta p) = G(s) \frac{ds}{d\Delta p}. \qquad (18)$$

We found [59, 63, 64, 71], both in simulations and experiments that the burst size distribution $G(s)$ is given by the scaling relation

$$G(s) = s^{-\tau}H\left(\frac{s}{s^*}\right), \qquad (19)$$

where the cutoff size of the bursts $s^*$ is [63].

$$s^* = n\kappa\Delta p^*/\bar{v} = \left( \frac{n\kappa}{\bar{v}N(p_c)} \right)^{\frac{\nu D}{1+\nu D}}. \qquad (20)$$

where $N(p_c)$ is the value of the normalized capillary pressure threshold distribution taken at the critical percolation threshold pressure $p_c$, and $\nu = 4/3$ is the exponent describing the divergence of the correlation length in percolation [76]. Since the capillary pressure threshold distribution $N(p_t)$ is normalized, the width of the distribution $\sigma \approx 1/N(p_c)$. Therefore the cutoff size of the bursts will increase with $\sigma$ according to Eq. 20 as $s^* \propto \sigma^{0.78}$. Martys et al. [77] derived the analytical form

$$\tau = 1 + \frac{D_f - 1/\nu}{D}, \qquad (21)$$

where $D$ and $D_f$ are, respectively, the fractal dimensions of the growing cluster and its front (seen as the green line Figure 1). Using the literature values $D = 1.82$ [22, 57] and $D_f = 4/3$ [23, 27, 58], we obtain $\tau = 1.32$, consistent with our experimental data shown in Figure 6.

We therefore get from Eqs 17, 18, 19, and 20

$$\mathcal{F}(\Delta p) = \Delta p^{-\tau} H\left(\frac{\Delta p}{\Delta p^*}\right) \left(\frac{n\kappa}{\bar{v}}\right)^{(1-\tau)}, \qquad (22)$$
<figure><img src="image_11.png" /><figcaption>**FIGURE 8** | Typical evolution of the pressure signal during an experiment. In this specific case, the liquid was withdrawn from the bottom at a rate of 1 ml/h and the model was positioned horizontally, i.e., $\theta = 0°$. The thick horizontal line denotes a characteristic average pressure during the bursts $\langle p \rangle = 210$ Pa. The left and right insets (shaded respectively in red and green) show zoomed in sections of the signal, where we can see the characteristic Haines jumps. The signal looks statistically similar in both insets, thus justifying the assumptions made in **Eq. 24**.</figcaption></figure>

where Eqs 20, 17 give the cutoff pressure $\Delta p^*$

$$ \Delta p^* = \frac{1}{N(p_c)} \left( \frac{n\kappa}{N(p_c)\bar{v}} \right)^{\frac{-1}{1+\gamma D}}. \quad (23) $$

**Figure 7** shows the scaling function $H(\Delta p/\Delta p^*)$ for various values of $\kappa$ from invasion percolation simulations [16] which confirms (22) and (23).

We will in the following use $\mathcal{F}(\Delta p)$ to calculate the work. Assume that the viscosity is low such that the total time of the bursts is short, and can be neglected compared to the total time considered. The work can then be calculated by averaging **Eq. 15**

$$ W = M n \kappa \langle p \rangle \int_{0}^{\infty} \Delta p \mathcal{F}(\Delta p) d\Delta p, \quad (24) $$

where $M$ is the total number of bursts during the time considered. Here we have assumed that the average pressure in each burst $(p_2 + p_3)/2$ in **Eq. 15** is independent of the capillary pressure drop of that burst and therefore can be averaged independently of $\Delta p$. The value of that average is $\langle p \rangle$. This assumption can be readily verified experimentally if we look at the typical temporal signature of the capillary pressure in one experiment, see **Figure 8**. The pressure signal in this figure corresponds to an experiment performed with a liquid withdrawal rate of 1 ml/h on a horizontal model, i. e, the angle $\theta = 0°$ (see **Figure 1**). The data is shown after the transient regime in which the capillary pressure grows from 0 to the fluctuating signal around a characteristic pressure $\langle p \rangle$. Notice that the typical size of the fluctuations in the pressure do not seem to depend strongly on the instantaneous value of the pressure, thus confirming the assumption made in the derivation of **Eq. 24**. This can be better seen from the insets in the shaded red and green

regions in the plot. We observe the typical Haines jumps signature characteristic of slow drainage processes. The signal on the left inset is statistically very similar to that on the right, thus further justifying the assumptions in **Eq. 24**), i.e., that the capillary pressure drop $\Delta p$ in a given burst does not depend strongly on the average pressure in the burst at steady state. The thick horizontal line in the plot denotes a characteristic pressure $\langle p \rangle = 210$ Pa.

By inserting the expression for $\mathcal{F}(\Delta p)$ from **Eq. 22** into **Eq. 24**, we have

$$ \begin{aligned} W = M \langle p \rangle \bar{v} \left( \frac{n\kappa}{\bar{v}} \right)^{2-\tau} & \int_{0}^{\infty} \Delta p^{1-\tau} H \left( \frac{\Delta p}{\Delta p^*} \right) d\Delta p = \\ & M \langle p \rangle \bar{v} \left( \frac{n\kappa}{\bar{v}} \right)^{2-\tau} \Delta p^{*2-\tau} I, \end{aligned} \quad (25) $$

where $I$ is the integral

$$ I = \int_{0}^{\infty} y^{1-\tau} H(y) dy. \quad (26) $$

Then by using **Eq. 23** we obtain

$$ W = M \bar{v} \langle p \rangle I \left( \frac{n\kappa}{\bar{v}N(p_c)} \right)^{\frac{(2-\tau)\gamma D}{1+\gamma D}}. \quad (27) $$

Assume that we consider the time $t_w$ the front needs to move a length $w$ corresponding to the width of the model. Within this time the invading fluid has invaded an area $a^2S$ of $S$ pores. The total number of bursts considered $M$ can be calculated as $S/\langle s \rangle$, where $\langle s \rangle$ is the average size of the bursts

$$ \langle s \rangle = I s^{*2-\tau} = I \left( \frac{\kappa n}{N(p_c)\bar{v}} \right)^{\frac{(2-\tau)\gamma D}{1+\gamma D}}. \quad (28) $$

where $I$ is the integral in (26). We then get the following simple expression for the work

$$ W = S \bar{v} \langle p \rangle. \quad (29) $$

This expression is the well-known relation of the work expressed as a pressure times a volume change $\langle p \rangle \Delta V$, where $\Delta V = S\bar{v} = Sa^2b = qt_w$ is the volume injected by the pump during the time $t_w$.

We can then calculate the total dissipation within $t_w$ by using **Equations 16, 29**

$$ \Phi = W - E_s = Sa^2b \langle p \rangle - \gamma_2 a^2 2S - C\gamma_2 ab 2S^{(D-1)/D} - \gamma_1 ab S_c, \quad (30) $$

where $C$ is a constant. Here the second term corresponds to the contribution to the surface energy from the top and the bottom interfaces (see discussion after **Eq. 15**) and the third term to the contribution from the side walls. The $D-1$ term in the exponent is due to the cut between the fractal invasion structure and the sidewall using one of Mandelbrot's rules of thumb [48, 49]. The last term corresponds to the surface energy connected to the contours of the trapped clusters $S_c$ seen in red in **Figure 4**. For large $S$ the third term can be neglected and

$$ \Phi \approx Sa^2b \langle p \rangle - \gamma_2 a^2 2S - \gamma_1 ab S_c. \quad (31) $$
<figure><img src="image_12.png" /><figcaption>FIGURE 9 | Energy balance showing how the width of the model w/a (horizontal axis) influences the applied work (thick black line), and how this work is split into energy dissipation (blue), surface energy associated to the top and bottom surfaces (green) and surface energy associated to the lateral surfaces (yellow).</figcaption></figure>

The length of the front $n$, seen as the green line in **Figure 4**, is given by

$$n \propto (w/a)^{D_f} \propto S^{D_f/D} \propto S_e^{D_f/D_e}, \quad (32)$$

where $D_f$ is the fractal dimension of the front (not included trapping) and $D_e$ is the fractal dimension of contours of the trapped clusters. From **Figure 4** we see that the measured fractal dimensions for $D_e$ and $D$ seem very close, effectively indistinguishable from one another considered the error bars. This result implies that $S \propto S_e$. In the experiment shown in **Figure 4** we have measured the proportionality relation to be $S_e = 1.32S$. We have also measured the prefactor linking $S$ and the window size $w/a$ (see relation (32)), obtaining the expression $S = 0.72 (w/a)^D$. By plugging these two equations into Eq. 31 for the dissipation $\phi$, we gain access to the full energy balance of the system, the resulting plot is shown in **Figure 9**. In this figure we see how the total applied work $W$ (thick black line) is split into the dissipated energy $\phi$ (blue) and the surface energy associated to creation of the top and bottom interfaces $E_{st} = \gamma_2 a^2 2S$ (green) and the lateral interface $E_{sl} = \gamma_1 a b S_e$ (yellow). The dissipation is calculated directly from Eq. 31. The total applied work, the surface energy associated to the top and bottom interfaces and that associated to the lateral interfaces are given respectively by the first, second and third terms in Eq. 31. Because all terms in (31) and 29 are proportional to $S$ the ratios $\phi/W$, $E_{sl}/W$ and $E_{st}/W$ are independent of $w/a$ for large systems.

In the computation of the energy balance in **Figure 9** using Eqs 31, 32, we used $\epsilon = 0.71$ measured from the experimental images (see **Figure 5**), $\gamma = 0.064N/m$, $\psi = 70^\circ$, $\gamma_1 = (1 - \epsilon)\gamma + \epsilon\gamma \cos(\psi) = 0.034N/m$, $\langle p \rangle = 210Pa$ and $D = 1.87$, see **Figure 4**. As previously stated, our assumptions for $\gamma_1$ and $\gamma_2$ correspond to the scenario in which no film of the defending wetting fluid is left behind covering the solid surfaces after the invasion by the nonwetting phase.

If instead of the quasi-two-dimensional system considered thus far we had a three-dimensional long square channel with a

cross sectional area of $w^2$, Eq. 29 for the work would remain the same, while Eq. 30 for the dissipation would be modified to

$$\Phi = Sa^3 \langle p \rangle - \gamma_1 a^2 S_e - 4 C \gamma_2 a^2 S^{(D-1)/D}, \quad (33)$$

where $C$ is a constant. When the system gets large the last term can be neglected compared to the two first terms such that

$$\Phi \approx Sa^3 \langle p \rangle - \gamma_1 a^2 S_e. \quad (34)$$

If the invaded surface area $S_e a^2$ times the characteristic pore size $a$ is proportional to the invaded volume $Sa^3$, the ratio $\phi/W$ and $E_{sl}/W$ will be independent of $w/a$.

## 4 STABILIZING FIELDS AND CROSSOVER LENGTHS

Let's now take into account a gravitational field and consider an experiment where a low density fluid with density $\rho_1$ and viscosity $\mu_1 = 0$ invades another fluid with a higher density $\rho_2$ and viscosity $\mu_2$ from above (See **Figure 1**). The mapping between the occupation probability in percolation theory $f$ and the capillary pressure $p$ is given by [76, 78].

$$f = \int_{0}^{p} N(p_t) dp_t. \quad (35)$$

We therefore have

$$f - f_c = \int_{p_c}^{p} N(p_t) dp_t, \quad (36)$$

where the critical occupation probability is $f_c$ and the critical percolation pressure is $p_c$. As before, $N(p_t)$ is the normalized capillary pressure threshold distribution. Now, by Taylor expanding $N(p_t)$ around $p_c$ to the lowest order in $p - p_c$ in Eq. 36 we get

$$f - f_c = N(p_c)(p - p_c). \quad (37)$$

It is reasonable to assume that the viscous pressure drop in the displaced fluid will depend linearly on the length scale since there are no trapped invading fluid clusters in the displaced fluid. Let us consider the capillary pressure at a height $x$. Since the gravitational field also depends linearly on the length scale $x$ we have

$$f - f_c = N(p_c) \left( \Delta \rho g a - \frac{q \mu_2 a}{k A} \right) (x - x_0)/a = F(x - x_0)/a, \quad (38)$$

where $x_0$ corresponds to a height along the front where the capillary pressure is at the percolation threshold. Here $k$ is the permeability felt by the displaced fluid, $A$ the cross-section area of the porous medium, and $\Delta \rho = \rho_2 - \rho_1$ the density difference between the fluids. The fluctuation number $F$

$$F = N(p_c) \left( \Delta \rho g a - \frac{q \mu_2 a}{k A} \right), \quad (39)$$

is a dimensionless number [61] which characterizes the gravitational and viscous fields and the capillary pressure
<figure><img src="image_13.png" /><figcaption>**FIGURE 10** | The dependence of the front width $\eta/a$ on the generalized Bond number $Fa/(\gamma N(p_c)) = Bo - Ca$. The blue dots correspond to data taken from [61] where the Bond number $Bo = \Delta\rho g a^2/\gamma = 0.154$ is kept fixed but the capillary number $Ca = q\mu_2 a^2/(\gamma k A)$ is changed by varying the flow rate $q$.</figcaption></figure>

fluctuations. Note that the width of the capillary pressure threshold distribution is $\sigma \approx 1/N(p_c)$, as the distribution is normalized. We will then use an assumption, first introduced by Sapoval [27, 57, 79, 80], that the correlation length in percolation will scale as the width of the front $\eta = (x - x_0) \propto \xi$, where $\xi$ is the percolation correlation length. This is based on the observation that the largest trapped cluster of wetting fluid is limited by the front width. From percolation theory, the correlation length $\xi$ is given by [76].

$$ (\xi/a) \propto (f - f_c)^{-\gamma}, \quad (40) $$

Using Sapoval's argument, and inserting Eq. 40 into Eq. 38 gives [26, 61].

$$ \eta/a \propto F^{\frac{\gamma}{1+\gamma}}, \quad (41) $$

Therefore, when $F > 0$, the front will be stabilized with a characteristic length scale $\eta$ [26, 27, 61]. This scaling behaviour is shown in **Figure 10**. The blue dots show experimental data taken from [61] where the authors analyzed the scaling of the invasion front in drainage by keeping the Bond number $Bo = \Delta\rho g a^2/\gamma = 0.154$ fixed but changing the capillary number $Ca$ by varying the flow rate $q$ (see **Figure 10** in Ref. [61] and experiments 1, 2, 3, 4, 5, 7, and 9 in Table I in the same reference).

Let us first assume that the gravitational effect is large enough such that it is the characteristic length scale $\eta$ that sets the width of the front $\eta < w$. Consider again a time scale $t_w$ corresponding to the time the front needs to move a length scale $w$. Dividing space into cells of size $\eta^2$, below which the fluid distribution is fractal, the number of pores $S$ and $S_e$ are given by

$$ S \propto \left(\frac{w}{\eta}\right)^2 \left(\frac{\eta}{a}\right)^D, \quad S_e \propto \left(\frac{w}{\eta}\right)^2 \left(\frac{\eta}{a}\right)^{D_e}, \quad (42) $$

and the work can be written as

$$ W = -g(\rho_2 h + \rho_1 (L-h))qt_w + \frac{q\mu_2 h}{kA}Sa^2b + Sa^2b\langle p \rangle, \quad (43) $$

where $h$ is the average position of the front from the lower outlet. The first term in (43) is the negative work due to the hydrostatic pressure, the second term the work due to the viscous pressure in the displaced fluid, and the last term is the work needed to continuously build up the interface energy. The work must be equal to the dissipation $\Phi$ plus the surface energy $E_s$ minus the change in gravitational potential energy $\Delta U_p$ (corresponding to the first term in Eq. 43)

$$ W = \Phi + E_s - \Delta U_p. \quad (44) $$

We then get the dissipation

$$ \Phi = \frac{q\mu_2 h}{\kappa A} Sa^2 b + Sa^2 b \langle p \rangle - \gamma_2 a^2 2S - C\gamma_2 ab 2S^{(D-1)/D} - \gamma_1 ab S_e^{D_e}, \quad (45) $$

which for large $S$ can be approximated as

$$ \Phi = \frac{q\mu_2 h}{\kappa A} Sa^2 b + Sa^2 b \langle p \rangle - \gamma_2 a^2 2S - \gamma_1 ab S_e^{D_e}. \quad (46) $$

Note that in the limit of very slow processes, the flow rate $q$ becomes negligible so the first term in this Eqs 45, 46 can be ignored, thus recovering Eqs 30, 31. In the scenario in which $\eta > w$, the relations for $S$ and $S_e$ will change to

$$ S \propto \left(\frac{w}{a}\right)^D, \quad S_e \propto \left(\frac{w}{a}\right)^{D_e}, \quad (47) $$

but we will still have the same expressions for $W$ and $\phi$. Again if $D_e = D$, the ratios $\phi/W$, $E_{sl}/W$ and $E_{st}/W$ will be independent of $w/a$ for large systems.

## 5 SATURATION BEHIND THE INVASION FRONT

The upscaling problem of calculating the large scale saturation involves identifying the cross-over length scale where the fluid distribution is no longer fractal. Let $L$ be the length of the porous model and assume $L > w$. The final saturation behind the front of the invading fluid and its dependence on the pressure across the model has been studied in Ref. [81]. In this study, we considered the volume of the invaded fluid in boxes with size corresponding to the width of the invading front. On length scales below this size, the structure of the invading fluid is fractal, while on length scales larger than this size, the structure is homogeneous. For sufficiently large $F$, when $\eta$ is the characteristic length scale of the front, the saturation $S_{nw}^F$ of the nonwetting fluid behind the front becomes [81].

$$ S_{nw}^F \propto \frac{\left(\frac{Lw^{d-1}}{\eta^d}\right) \left(\frac{\eta}{a}\right)^D a^d}{Lw^{d-1}} = \left(\frac{\eta}{a}\right)^{D-d}, \quad (48) $$

where $d$ is the spatial dimension (2 or 3). Hence using Eq. 41

$$ S_{nw}^F \propto F^{-\gamma(D-d)/(\gamma+1)}. \quad (49) $$

However, when the fluctuation number $F$ is sufficient small, $\eta > w$, and the width of the model $w$ will be the characteristic length scale in the problem. Then
<figure><img src="image_14.png" /><figcaption>**FIGURE 11** | The final saturation as a function of the generalized Bond number $Fa/(γN(p_c)) = Bo - Ca$ is plotted both for experimental (black dots) and numerical results (blue stars) produced via invasion percolation with a gravitational field, together with the predicted result Eq. 49 (red dash-dotted line) in two dimensions. The capillary number $Ca = qμ_2a^2/(γkA) = 1.2 \cdot 10^{-4}$ is kept fixed while the Bond number $Bo = Δρga^2/γ$ is changed in the experiments. Numerical and experimental data from [81].</figcaption></figure>

$$S_{nw}^F \propto \left(\frac{w}{a}\right)^{(D-d)} \qquad (50)$$

**Figure 11** shows a two-dimensional invasion percolation simulation with a gravitational field together with drainage experiments performed by Ayaz et al. [81]. The red dash-dotted line confirms the predictions of the theoretical scaling in Eq. 49. We make the important remark here that in the equations considered above we have neglected the boundary effects at the inlet and outlet since $L/w \gg 1$. These boundary effects are however important when the characteristic length scale of the front $w$ is of the same order as the length $L$ of the system [74].

In section IV we considered an invading fluid with a negligible viscosity $μ_1 = 0$. Let us now instead consider the scenario in which this fluid has a non-negligible viscosity $μ_1$. In that case the fluctuation number will be

$$F = N(p_c) \left( \Delta \rho g a + \frac{q \mu_1 a}{k_1 A} - \frac{q \mu_2 a}{k A} \right), \qquad (51)$$

where $k_1$ is an effective permeability that depends on the characteristic crossover length scale of the system which is either $η$ or $w$. The expression (51) is meaningful only if the pressure drop in the invading fluid depends linearly on the length scale along the average flow direction. Network simulations in two dimensions however, show a linear change in the capillary pressure due to viscous flow [59, 60]. This is due to the loop-less strands in the invading nonwetting fluid. It is of great interest to know if this is also valid in three dimensions.

The structure left behind the front is a fractal structure with clusters of all length scales up to the characteristic length scale $η$ or $w$ [27, 58]. Therefore if $η < w$, the permeability $k_1(F)$ will be a

function of $F$. If we know this functional dependence we can solve the equation

$$F = N(p_c) \left( \Delta \rho g a + \frac{q \mu_1 a}{k_1 (F) A} - \frac{q \mu_2 a}{k A} \right), \qquad (52)$$

with respect to $F$ and thereby find the characteristic length scale $η$ from

$$η/a \propto F^{\frac{2}{1+\gamma}}. \qquad (53)$$

We then also know the saturation through Eq. 48. A goal for future experiments/simulations/theory should therefore be to measure/predict the permeability dependence of the invading fluid as function of the characteristic length scale $η$.

# 6 CONCLUSION

In this paper, we discussed the importance of capillary fluctuations in porous media, as well as the characteristic length scales set by the competition between capillary fluctuations and external fields, gravitational or viscous. We focused on fluid fronts that are stabilized by gravitational and viscous fields. The fluctuation number $F$, which describes the scaling of the front width $η$, was introduced. When considering viscous and gravitational fields, 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 experiments. 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 $η$ is thus of primary importance in defining a relevant Representative Elementary Volume (REV) for an average Darcy description of the two-phase flow problem.

We also discussed the energy dissipation caused by out-of-equilibrium Haines jumps and calculated the total elastic energy released by those jumps using the pressure drops connected to the jumps' known scaling behavior. When the invasion front is in a steady state regime, we find that the elastic energy released by the jumps corresponds to the work applied to the system, as expected. However, this description also explains the local activity and dissipation caused by local Haines jumps, as well as how it is related to dissipation on a larger scale.

We calculated the energy dissipation and saturation left behind the invasion front in addition to the work done by the pump to push the fluid front forward. These quantities are affected by the generalized fluctuation number, as well as the fractal dimension of the invasion structure, the fractal dimension of the contour of the trapped clusters, and the exponent $ν$, which describes the divergence of the correlation length in percolation as it approaches the critical percolation threshold $f_c$. The work $W$, dissipation $Φ$, and surface energy $E_s$ were discovered to be proportional to the number of pores invaded $S$, with fractal scaling on length scales smaller than the characteristic length scale, $η$ or $w$, but scaling with the spatial dimension on larger
length scales. Hence, the ratio of these energies stays constant as the system size increases.

The theory described here is a bottom-up approach that allows the problem to be scaled up from small to large scales. The theoretical results were compared to and found to be consistent with quasi-two-dimensional experiments. Because of the small temperature increase involved in heat dissipation from the Haines jumps (mK or less), it is very challenging to measure these temperature fluctuations. Such a measurement remains to be done. An experimental verification of the theoretical scaling of the invasion front and the saturation behind the front for three-dimensional porous media is also of great interest for future work [81]. Experiments involving changing the gravitational effects in a geophysical centrifuge could be one way to accomplish this. The porous medium in a real rock is typically not homogeneous, in contrast to the porous media considered here. As a result, it is of great interest to extend the theory to include also the case of inhomogeneous porous media, for example with a porosity gradient.

## DATA AVAILABILITY STATEMENT

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## REFERENCES

1. Anderson R, Zhang L, Ding Y, Blanco M, Bi X, Wilkinson DP. A Critical Review of Two-phase Flow in Gas Flow Channels of Proton Exchange Membrane Fuel Cells. *J Power Sourc* (2010) 195:4531-53. doi:10.1016/j.jpowsour.2009.12.123

2. Hwang SH, Kim C, Song H, Son S, Jang J. Designed Architecture of Multiscale Porous Tio2 Nanofibers for Dye-Sensitized Solar Cells Photoanode. *ACS Appl Mater Inter* (2012) 4:5287-92. doi:10.1021/am301245s

3. Cantwell WJ, Morton J. The Impact Resistance of Composite Materials - a Review. *Composites* (1991) 22:347-62. doi:10.1016/0010-4361(91)90549-v

4. Gibson P, Rivin D, Kendrick C, Schreuder-Gibson H. Humidity-dependent Air Permeability of Textile Materials1. *Textile Res J* (1999) 69:311-7. doi:10.1177/004051759906900501

5. Lake L. *Enhanced Oil Recovery*. Englewood Cliffs, N.J.: Prentice-Hall (1989).

6. Shabani Afrapoli M, Alipour S, Torsaeter O. Fundamental Study of Pore Scale Mechanisms in Microbial Improved Oil Recovery Processes. *Transp Porous Med* (2011) 90:949-64. doi:10.1007/s11242-011-9825-7

7. Yan J, Luo X, Wang W, Toussaint R, Schmittbuhl J, Vasseur G, et al. An Experimental Study of Secondary Oil Migration in a Three-Dimensional Tilted Porous Medium. *Bulletin* (2012) 96:773-88. doi:10.1306/09091110140

8. Vasseur G, Luo X, Yan J, Loggia D, Toussaint R, Schmittbuhl J. Flow Regime Associated with Vertical Secondary Migration. *Mar Pet Geology* (2013) 45: 150-8. doi:10.1016/j.marpetgeo.2013.04.020

9. Guymon G. *Unsaturated Zone Hydrology*. Englewood Cliffs, N.J.: Prentice-Hall (1994).

10. Bear J. *Dynamics of Fluids in Porous media*. New York: Elsevier (1972).

11. Bear J, Verruijt A. *Modeling Groundwater Flow and Pollution*. Dordrecht Boston Norwell, MA, U.S.A.: D. Reidel Pub. Co. Sold and distributed in the U.S.A. and Canada by Kluwer Academic Publishers (1987).

12. Jellali S, Muntzer P, Razakarisoa O, Schäfer G. Large Scale experiment on Transport of Trichloroethylene in a Controlled Aquifer. *Transp Porous Media* (2001) 44:145-63. doi:10.1023/A:1010652230922

## AUTHOR CONTRIBUTIONS

KM initiated the work, did the main part of the theory, and wrote the first draft of the article. MM and KM designed the experiment. MM performed the experiment and the data analyzes. MM, EF, RT, AH and KM were involved in the discussions that led to the theoretical analysis presented. All authors read critically and participated to the writing of the manuscript.

## FUNDING

This work was partly supported by the Research Council of Norway through its Centers of Excellence funding scheme, project number 262 644.

## ACKNOWLEDGMENTS

We thank Einar Hinrichsen, Jens Feder, Torstein Jøssang, Amnon Aharony, Paul Meakin, Jean Schmittbuhl, Gerhard Schäfer, Liv Furuberg, Alexander Birovljev, Yves Meheust, Per Arne Rikvold, Harold Auradou, Olav Inge Frette, Vidar Frette, Grunde Løvoll, Bjørnar Sandnes, Mihailo Jankov, Fredrik Kvalheim Eriksen, Monem Ayaz, and Tom Vincent-Dospital for fruitful discussions.

13. Nsir K, Schäfer G, Roupert Rd. C, Razakarisoa O, Toussaint R. Laboratory Experiments on DNAPL Gravity Fingering in Water-Saturated Porous media. *Int J Multiphase Flow* (2012) 40:83-92. doi:10.1016/j.ijmultiphaseflow.2011.12.003

14. Gagné J. *The Physics of Filter Coffee*. Scott Rao Coffee Books (2021).

15. Chen J-D, Wilkinson D. Pore-scale Viscous Fingering in Porous media. *Phys Rev Lett* (1985) 55:1892-5. doi:10.1103/physrevlett.55.1892

16. Måløy KJ, Feder J, Jøssang T. Viscous Fingering Fractals in Porous media. *Phys Rev Lett* (1985) 55:2688-91.

17. Weitz DA, Stokes JP, Ball RC, Kushnick AP. Dynamic Capillary Pressure in Porous media: Origin of the Viscous-Fingering Length Scale. *Phys Rev Lett* (1987) 59:2967-70. doi:10.1103/physrevlett.59.2967

18. Lenormand R, Touboul E, Zarcone C. Numerical Models and Experiments on Immiscible Displacements in Porous media. *J Fluid Mech* (1988) 189:165-87. doi:10.1017/s0022112088000953

19. Løvoll G, Mèheust Y, Toussaint R, Schmittbuhl J, Måløy KJ. Growth Activity during Fingering in a Porous Hele Shaw Cell. *Phys Rev E* (2004) 70:026301.

20. Tallakstad KT, Knudsen HA, Ramstad T, Løvoll G, Måløy KJ, Toussaint R, et al. Steady-state Two-phase Flow in Porous media: Statistics and Transport Properties. *Phys Rev Lett* (2009) 102:074502. doi:10.1103/PhysRevLett.102.074502

21. Tallakstad KT, Løvoll G, Knudsen HA, Ramstad T, Flekkøy EG, Måløy KJ. Steady-state, Simultaneous Two-phase Flow in Porous media: An Experimental Study. *Phys Rev E Stat Nonlin Soft Matter Phys* (2009) 80: 036308. doi:10.1103/PhysRevE.80.036308

22. Lenormand R, Zarcone C. Invasion Percolation in an Etched Network: Measurement of a Fractal Dimension. *Phys Rev Lett* (1985) 54:2226-9. doi:10.1103/physrevlett.54.2226

23. Furuberg L, Feder J, Aharony A, Jøssang T. Dynamics of Invasion Percolation. *Phys Rev Lett* (1988) 61:2117. doi:10.1103/physrevlett.61.2117

24. Moura M, Måløy KJ, Flekkøy EG, Toussaint R. Verification of a Dynamic Scaling for the Pair Correlation Function during the Slow Drainage of a Porous Medium. *Phys Rev Lett* (2017) 119:154503. doi:10.1103/physrevlett.119.154503

25. Xiao B, Huang Q, Chen H, Chen X, Long G. A Fractal Model for Capillary Flow through a Single Tortuous Capillary with Roughened Surfaces in
Fibrous Porous media. Fractals (2021) 29:2150017. doi:10.1142/S0218348X21500171

26. Wilkinson D. Percolation Model of Immiscible Displacement Displacement in the Presence of Buoyancy Forces. *Phys Rev A* (1984) 30:520. doi:10.1103/physreva.30.520

27. Birovljev A, Furuberg L, Feder J, Jøssang T, Måløy KJ, Aharony A. Gravity Invasion Percolation in Two Dimensions: Experiment and Simulation. *Phys Rev Lett* (1991) 67:584. doi:10.1103/physrevlett.67.584

28. Frette V, Feder J, Jøssang T, Meakin P. Buoyancy Fluid Migration in Porous media. *Phys Rev Lett* (1992) 68:3164. doi:10.1103/physrevlett.68.3164

29. Wagner G, Birovljev A, Meakin P, Feder J, Jøssang T. Fragmentation and Migration of Invasion Percolation Clusters:experiments and Simulations. *Phys Rev E* (1997) 55:7015. doi:10.1103/physreve.55.7015

30. Auradou H, Måløy KJ, Schmittbuhl J, Hansen A, Bideau D. Competition between Correlated Buoyancy and Uncorrelated Capillary Effects during Drainage. *Phys Rev E* (1999) 60:7224. doi:10.1103/physreve.60.7224

31. Muharrik M, Suekane T, Patmonoaji A. Effect of Buoyancy on Fingering Growth Activity in Immiscible Two-phase Flow Displacement. *J Fluid Sci Techn* (2018) 13:17. doi:10.1299/jfst.2018jfst0006

32. Zhao B, Mac Minn CW, Juanes R. Wettability Control on Multiphase Flow in Patterned Microfluidics. *PNAS* (2016) 113:10251. doi:10.1073/pnas.1603387113

33. Cottin C, Bodiguel H, Colin A. Influence of Wetting Conditions on Drainage in Porous media: A Microfluidic Study. *Phys Rev E* (2011) 84:026311. doi:10.1103/PhysRevE.84.026311

34. Holtzman R, Segre E. Wettability Stabilizes Fluid Invasion into Porous media via Nonlocal, Cooperative Pore Filling. *Phys Rev Lett* (2015) 115:164501. doi:10.1103/physrevlett.115.164501

35. Holtzman R, Dentz M, Planet R, Ortin J. The Origin of Hysteresis and Memory of Two-phase Flow in Disordered media. *Commun Phys* (2020) 3:222. doi:10.1038/s42005-020-00492-1

36. de Gennes PG, Brochard-Wyart F, Quéré D. Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves. New York: Springer (2004).

37. Primkulov BK, Talman S, Khaleghi K, Rangriz Shokri A, Chalaturnyk R, Zhao B, et al. Quasistatic Fluid-Fluid Displacement in Porous media: Invasion-Percolation through a Wetting Transition. *Phys Rev Fluids* (2018) 3:104001. doi:10.1103/physrevfluids.3.104001

38. Primkulov BK, Pahlavan AA, Fu X, Zhao B, MacMinn CW, Juanes R. Signatures of Fluid-Fluid Displacement in Porous media: Wettability, Patterns and Pressures. *J Fluid Mech* (2019) 875:R4. doi:10.1017/jfm.2019.554

39. Rabbani HS, Or D, Liu Y, Lai C-Y, Lu NB, Datta SS, et al. Suppressing Vicious Fingering in Structures Porous media. *PNAS* (2018) 115:4833. doi:10.1073/pnas.1800729115

40. Lu NB, Browne CA, Amchin DB. Controlling Capillary Fingering Using Pore Size Gradients in Disordered media. *Phys Rev Fluids* (2019) 4:084303. doi:10.1103/physrevfluids.4.084303

41. Lenormand R. Flow through Porous media: Limits of Fractal Pattern. *Proc R Soc Lond A* (1989) 423:159.

42. Payatakes AC, Ng KM, Flumerfelt RW. Oil Ganglion Dynamics during Immiscible Displacement: Model Formulation. *AIChE J* (1980) 26:430–43. doi:10.1002/aic.690260315

43. Sandnes B, Knudsen HA, Måløy KJ, Flekkøy EG. Labyrinth Patterns in Confined Granular-Fluid Systems. *Phys Rev Lett* (2007) 99:038001. doi:10.1103/PhysRevLett.99.038001

44. Sandnes B, Flekkøy EG, Måløy KJ, See H. Patterns and Flow in Frictional Fluid Dynamics. *Nat Commun* (2011) 2:288. doi:10.1038/ncomms1289

45. Odier C, Levaché B, Santanach-Carreras E, Bartolo D. Forced Imbibition in Porous media: A Fourfold Scenario. *Phys Rev Lett* (2017) 119:208005. doi:10.1103/physrevlett.119.208005

46. Ferer M, Sams WN, Geisbrecht RA, Smith DH. Crossover from Fractal to Compact Growth from Simulations of Two-phase Flow with Finite Viscosity Ratio in Two-Dimensional Porous media. *Phys Rev E* (1993) 47:2713. doi:10.1103/physreve.47.2713

47. Ferer M, Bromhal GS, Smith DH. Pore-level Modeling of Drainage: Crossover from Invasion Percolation to Compact Flow. *Phys Rev E* (2003) 67:051601. doi:10.1103/PhysRevE.67.051601

48. B Mandelbrot B. *The Fractal Geometry of Nature*. San Francisco: W. H. Freeman (1982).

49. Feder J. Fractals. New York: Plenum (1988).

50. Yu B. Analysis of Flow in Fractal Porous media. *Appl Mech Rev* (2008) 61: 050801. doi:10.1115/1.2955849

51. Xiao B, Wang W, Zhang X, Long G, Fan J, Chen H, et al. A Novel Fractal Solution for Permeability and Kozeny-Carman Constant of Fibrous Porous media Made up of Solid Particles and Porous Fibers. *Powder Techn* (2019) 349: 92. doi:10.1016/j.powtec.2019.03.028

52. Liang M, Fu C, Xiao B, Luo L, Wang Z. A Fractal Study for the Effective Electrolyte Diffusion through Charged Porous media. *Int J Heat Mass Transfer* (2019) 137:365. doi:10.1016/j.ijheatmasstransfer.2019.03.141

53. Armand M, Axmann P, Bresser D, Copley M, Edström K, Ekberg C, et al. Lithium-ion Batteries – Current State of the Art and Anticipated Developments. *J Power Source* (2020) 479:228708. doi:10.1016/j.jpowsour.2020.228708

54. Darcy H. *Les fontaines publiques de dijon*. Paris: Dalmont (1856).

55. de Gennes PG, Guyon E. Lois Generales Pour Linjections Dún Fluide Dans Un Milieu Poreux Aleatoire. *J Mech* (1978) 17:403.

56. Chandler R, Koplik J, Lerman K, Willemsen D. Capillary Displacement and Percolation in Porous media. *J Fluid Mech* (1982) 119:249. doi:10.1017/s0022112082001335

57. Wilkinson D, Willemsen JF. Invasion Percolation: a New Form of Percolation Theory. *J Phys A: Math Gen* (1983) 16:3365. doi:10.1088/0305-4470/16/14/028

58. Frette OI, Måløy KJ, Schmittbuhl J. Immiscible Displacement of Viscosity-Matched Fluids in Two-Dimensional Porous media. *Phys Rev E* (1997) 57: 2969. doi:10.1103/physreve.55.2969

59. Aker E, Jørgen Måløy K, Hansen A. Viscous Stabilization of 2D Drainage Displacements with Trapping. *Phys Rev Lett* (2000) 84(20):4589–92. doi:10.1103/PhysRevLett.84.4589

60. Aker E, Måløy KJ, Hansen A. Dynamics of Stable Viscous Displacement in Porous media. *Phys Rev E* (2000) 61:2936. doi:10.1103/physreve.61.2936

61. Meheyst Y, Løvoll G, Måløy KJ, Schmittbuhl J. Interface Scaling in a Two-Dimensional Porous Medium under Combined Viscous, Gravity, and Capillary Effects. *Phys Rev E* (2002) 66:051603. doi:10.1103/PhysRevE.66.051603

62. Haines WB. Studies in the Physical Properties of Soil. The Hysteresis Effect in Capillary Properties, and the Modes of Moisture Distribution Associated Therewith. *J Agric Sci* (1930) 20:97. doi:10.1017/s002185960008864x

63. Måløy KJ, Furuberg L, Feder J, Jøssang T. Dynamics of Slow Drainage in Porous media. *Phys Rev Lett* (1992) 68:2161.

64. Furuberg L, Måløy KJ, Feder J. Intermittent Behaviour in Slow Drainage. *Phys Rev E* (1996) 53:966. doi:10.1103/physreve.53.966

65. Moebius F, Or D. Pore Scale Dynamics Underlying the Motion of Drainage Fronts in Porous media. *Water Resour Res* (2014) 50:8441. doi:10.1002/2014wr015916

66. Moebius F, Or D. Inertial Forces Affect Fluid Front Displacement Dynamics in a Pore-Throat Network Model. *Phys Rev E* (2014) 90:023019. doi:10.1103/PhysRevE.90.023019

67. Moebius F, Canone D, Or D. Characteristics of Acoustic Emissions Induced by Fluid Front Displacement in Porous media. *Water Resour Res* (2012) 48: W11507. doi:10.1029/2012wr012525

68. Berg S, Ott H, Klapp SA, Schwing A, Neiteler R, Brussee N, et al. Real-time 3d Imaging of haines Jumps in Porous media Flow. *PNAS* (2013) 110:3755–9. doi:10.1073/pnas.1221373110

69. Bultreys T, Boone MA, Boone MN, De Schryver T, Masschaele B, Van Loo D, et al. Real-time Visualization of haines Jumps in sandstone with Laboratory-Based Microcomputed Tomography. *Water Resour Res* (2015) 51:8668. doi:10.1002/2015wr017502

70. Zacharoudiou I, Boek ES, Crawshaw J. The Impact of Drainage Displacement Patterns and haines Jumps on Co2 Storage Efficiency. *Scientific Rep* (2018) 8: 15561. doi:10.1038/s41598-018-33502-y

71. Moura M, Måløy KJ, Toussaint R. Critical Behaviour in Porous media Flow. *Europhys Lett* (2017) 118:14004. doi:10.1209/0295-5075/118/14004

72. Berg CF, Slotte PA, Khanamiri HH. Geometrical Derived Efficiency of Slow Immiscible Displacement in Porous media. *Phys Rev E* (2020) 102:033113. doi:10.1103/PhysRevE.102.033113
73. Hele-Shaw HS. The Flow of Water. *Nature* (1898) 58:34. doi:10.1038/058034a0

74. Moura M, Florentino EA, Måløy KJ, Schäfer G, Toussaint R. Impact of Sample Geometry on the Measurement of Pressure-Saturation Curves: Experiments and Simulations. *Water Resour Res* (2015) 51. doi:10.1002/2015WR017196

75. Hinrichsen EL, Feder J, Jøssang T. Geometry of Random Sequential Adsorption. *Jornal Stat Phys* (1986) 44:793. doi:10.1007/bf01011908

76. Stauffer D. *Introduction to Percolation Theory*. London Bristol, PA: Taylor & Francis (1994).

77. Martys N, Robbins MO, Cieplak M. Scaling Relations for Interface Motion through Disordered media: Application to Two-Dimensional Fluid Invasion. *Phys Rev B* (1991) 44:12294. doi:10.1103/physrevb.44.12294

78. Roux S, Guyon E. Temporal Development of Invasion Percolation. *J Phys A* (1989) 22:3693. doi:10.1088/0305-4470/22/17/034

79. Sapoval E, Rosso M, Gouyet JF. The Fractal Nature of a Diffusion Front and the Relation to Percolation. *J Phys (France) Lett* (1985) 46:L149. doi:10.1051/jphyslet:01985004604014900

80. Gouyet JF, Rosso M, Sapoval B. Fractal Structure of Diffusion and Invasion Fronts in Three-Dimensional Lattices through the Gradient Percolation Approach. *Phys Rev B* (1988) 37:1832. doi:10.1103/physrevb.37.1832

81. Ayaz M, Toussaint R, Schäfer G, Måløy KJ. Gravitational and Finite-Size Effects on Pressure Saturation Curves during Drainage. *Water Resour Res* (2020) 56:e2019WR026279. doi:10.1029/2019wr026279

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Copyright © 2021 Måløy, Moura, Hansen, Flekkøy and Toussaint. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.