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

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.

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 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 κ 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 = ⟨p⟩ΔV, where ⟨p⟩ is the average pressure across the model and Δ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 Φ, 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 η, the spatial scaling of S is fractal and for larger length scales the scaling becomes proportional to the system size instead. We then show that we can explicitly calculate W, E_s, and Φ 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 Φ/W are both independent of system size.

FIGURE 1

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₀ sin(θ) along the x direction of the model, where g₀ = 9.82 m/s². The average position of the front is h and the width of the front is 2η. 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.

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 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η.

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 ν = 4.25 ⋅ 10^–5 m²/s, ρ = 1.205 g/cm³ and γ = 0.064 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 ψ = 70°. 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₀ to one with a radius of curvature r = r₀—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₀ to the front defined by the radius of curvature r = r₀—dr. The longest of these distances we denote as dx_i. The distance dx_i to lowest order in dr is given by

$d{x}_i={\left.\frac{d{x}_i}{dr}\right|}_{r_0}dr={\alpha }_{i}dr,(1)$

where

${\alpha }_i(\psi ,p_0)={\left.\frac{d{x}_i}{dr}\right|}_{r_0}(2)$

depends on the wetting angle ψ (See Figure 2) and the capillary pressure p₀ given by the Young-Laplace equation p = γ(1/r + 1/r₁) [10] for r = r₀. Here r is the in-plane, and r₁ the out of plane radius of curvature assumed to be constant in the quasi-two-dimensional experiments.

FIGURE 2

FIGURE 2. The figure illustrates a nonwetting fluid invading a wetting fluid at a single pore-throat at different radius of curvature r₀ and r = r₀ − dr. The distance moved between the two front positions is dx_i and ψ is the wetting angle.

We can then calculate the corresponding increase in capillary pressure

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

Then by using Eqs 1, 3 we get

$dp=\frac{\gamma }{r_o^2{\alpha }_i}d{x}_i.(4)$

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

$d{v}_i={\mathcal{A}}_i(\psi ,p_0)d{x}_i,(5)$

where ${\mathcal{A}}_i(\psi ,p_0)$ is the surface area in pore i which is function of the wetting angle and the capillary pressure p₀. We then get

$dp=\frac{\gamma }{r_o^2{\mathcal{A}}_i{\alpha }_i}d{v}_i,(6)$

such that

$d{v}_i={\kappa }_{i}dp,(7)$

where the capacitive volume

${\kappa }_i(\psi ,p_0)=\frac{r_0^2{\mathcal{A}}_i(\psi ,p_0){\alpha }_i(\psi ,p_0)}{\gamma },(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{\mathcal{A}}_{i}d{x}_i,(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}^{\ast }d{p}^{\ast }=U_i(0)+\frac{{\kappa }_i}{2}{p}^2.(10)$

We have assumed p₀ = 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,(11)$

where

$\kappa =\left(\sum _{i=1}^n{\kappa }_i\right)/n,(12)$

where κ_i is the capacitive volume of pore i so that κ is the average capacitive volume over the n interface pores. The capillary pressure p₂ at the time t₂ and the capillary pressure p₁ at the time t₁ are related by

$p_2-p_1=\frac{q(t_2-t_1)}{n\kappa },(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₁ at t₁ to p₂ at t₂ 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₁) to U_t (t₂)

$\begin{array}{l}\hfill U_t(t_2)-U_t(t_1)=W=n\kappa {\int }_{p_1}^{p_2}{p}^{\ast }d{p}^{*}=\\ \hfill \frac{n\kappa }{2}(p_2^2-p_1^2).\end{array}(14)$

FIGURE 3

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.

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₂ at a capillary pressure p₂ = p_t and that, after the burst, the system reaches another equilibrium state at capillary pressure p₃. 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

$\begin{array}{l}\hfill U_t(p_2)-U_t(p_3)=\frac{\kappa n(t_2)}{2}(p_2^2-p_3^2)=\\ \hfill \kappa n(t_2)(p_2-p_3)(p_2+p_3)/2,\end{array}(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 Δp = p₂ − p₃ is the pressure drop in the burst and (p₂ + p₃)/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γ₂a²s and one from the interface contour aγ₁bs_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 as_e is the length of the invasion contour of the burst (including trapped clusters). The factor 2 in the term 2γ₂a²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, γ₁ and γ₂, 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 ± 0.10, D_e = 1.84 ± 0.10 and D_f = 1.45 ± 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 δ < 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).

FIGURE 4

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_e of the invasion contours and the fractal dimension D_f of the main invasion front (green line). The dotted and dashed lines have slopes of −2 and −1 respectively.

When a fluid film is left behind the invasion front the surface tension γ₁ = γ₂ = γ which is the surface tension between the two fluids. However when no film is left behind the invasion front, γ₂ = γ_ns − γ_ws, where γ_ns is the interface tension between the nonwetting fluid and the solid and γ_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 γ₂ = γ cos(ψ), where ψ 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) 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 ϵ = 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, ϵ 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 γ₁ = (1 − ϵ)γ + ϵ(γ_ns − γ_ws) for the total internal interfacial contours. Using again the Young equation (10) we can rewrite this expression as γ₁ = (1 − ϵ)γ + ϵγ cos(ψ).

FIGURE 5

FIGURE 5. Tracking of the internal interfacial contours between the nonwetting and wetting phases (orange) and the nonwetting and solid phases (green).

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 ϕ plus the surface energy needed to create the new interface E_s, such that

$W=\mathrm{\Phi }+E_s.(16)$

It is important to note that W, Φ, and E_s are average quantities that do not account for fluctuations, and that ϕ 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}(\mathrm{\Delta }p)$ during the bursts since

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

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

$\mathcal{F}(\mathrm{\Delta }p)=G(s)\frac{ds}{d\mathrm{\Delta }p}.(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^{\ast }}\right),(19)$

where the cutoff size of the bursts s∗ is [63].

$s^{\ast }=n\kappa \mathrm{\Delta }{p}^{\ast }/\bar{v}={\left(\frac{n\kappa }{\bar{v}N(p_c)}\right)}^{\frac{\nu D}{1+\nu D}}.(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 ν = 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 σ ≈ 1/N (p_c). Therefore the cutoff size of the bursts will increase with σ according to Eq. 20 as s^∗ ∝ σ^0.78. Martys et al. [77] derived the analytical form

$\tau =1+\frac{D_f-1/\nu }{D},(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 τ = 1.32, consistent with our experimental data shown in Figure 6.

FIGURE 6

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].

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

$\mathcal{F}(\mathrm{\Delta }p)=\mathrm{\Delta }{p}^{-\tau }H\left(\frac{\mathrm{\Delta }p}{\mathrm{\Delta }{p}^{\ast }}\right){\left(\frac{n\kappa }{\bar{v}}\right)}^{(1-\tau )},(22)$

where Eqs 20, 17 give the cutoff pressure Δp∗

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

Figure 7 shows the scaling function H (Δp/Δp∗) for various values of κ from invasion percolation simulations [16] which confirms (22) and (23).

FIGURE 7

FIGURE 7. The dependence of the crossover function H (Δp/Δp*) of the pressure jump distribution in a modified invasion percolation simulation [63]. For each value of κ, the points represent averages over five independent simulations on 200 × 1,500 lattices. The scaling exponent 0.30 ≈ 1/(1 + νD). Data is taken from [63].

We will in the following use $\mathcal{F}(\mathrm{\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=Mn\kappa ⟨p⟩{\int }_0^{\infty }\mathrm{\Delta }p\mathcal{F}(\mathrm{\Delta }p)d\mathrm{\Delta }p,(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₂ + p₃)/2 in Eq. 15 is independent of the capillary pressure drop of that burst and therefore can be averaged independently of Δp. The value of that average is ⟨p⟩. 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 θ = 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 ⟨p⟩. 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 Δ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 ⟨p⟩ = 210 Pa.

FIGURE 8

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., θ = 0°. The thick horizontal line denotes a characteristic average pressure during the bursts ⟨p⟩ = 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.

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

$\begin{array}{l}\hfill W=M\langle p\rangle \bar{v}{\left(\frac{n\kappa }{\bar{v}}\right)}^{2-\tau }{\int }_0^{\infty }\mathrm{\Delta }{p}^{1-\tau }H\left(\frac{\mathrm{\Delta }p}{\mathrm{\Delta }{p}^{\ast }}\right)d\mathrm{\Delta }p=\\ \hfill M\langle p\rangle \bar{v}{\left(\frac{n\kappa }{\bar{v}}\right)}^{2-\tau }{\mathrm{\Delta }{p}^{\ast }}^{2-\tau }I,\end{array}(25)$

where I is the integral

$I={\int }_0^{\infty }{y}^{1-\tau }H(y)dy.(26)$

Then by using Eq. 23 we obtain

$W=M\bar{v}⟨p⟩I{\left(\frac{n\kappa }{\bar{v}N(p_c)}\right)}^{\frac{(2-\tau )\nu D}{1+\nu D}}.(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²S of S pores. The total number of bursts considered M can be calculated as S/⟨s⟩, where ⟨s⟩ is the average size of the bursts

$⟨s⟩=I{s^{\ast }}^{2-\tau }=I{\left(\frac{\kappa n}{N(p_c)\bar{v}}\right)}^{\frac{(2-\tau )\nu D}{1+\nu D}}.(28)$

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

$W=S\bar{v}⟨p⟩.(29)$

This expression is the well-known relation of the work expressed as a pressure times a volume change ⟨p⟩ΔV, where $\mathrm{\Delta }V=S\bar{v}=S{a}^{2}b=q{t}_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

$\begin{array}{ll}\hfill \mathrm{\Phi }=W-E_s=& \hfill S{a}^{2}b⟨p⟩-{\gamma }_2{a}^{2}2S-C{\gamma }_{2}ab2S^{(D-1)/D}-{\gamma }_{1}ab{S}_e,\end{array}(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_e seen in red in Figure 4. For large S the third term can be neglected and

$\mathrm{\Phi }\approx S{a}^{2}b⟨p⟩-{\gamma }_2{a}^{2}2S-{\gamma }_{1}ab{S}_e.(31)$

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},(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 ∝ 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 ϕ, 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 ϕ (blue) and the surface energy associated to creation of the top and bottom interfaces E_st = γ₂a²2S (green) and the lateral interface E_sl = γ₁abS_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 ϕ/W, E_sl/W and E_st/W are independent of w/a for large systems.

FIGURE 9

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).

In the computation of the energy balance in Figure 9 using Eqs 31, 32, we used ϵ = 0.71 measured from the experimental images (see Figure 5), γ = 0.064N/m, ψ = 70°, γ₁ = (1 − ϵ)γ + ϵγ cos(ψ) = 0.034N/m, ⟨p⟩ = 210Pa and D = 1.87, see Figure 4. As previously stated, our assumptions for γ₁ and γ₂ 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², Eq. 29 for the work would remain the same, while Eq. 30 for the dissipation would be modified to

$\mathrm{\Phi }=S{a}^3⟨p⟩-{\gamma }_1{a}^2{S}_e-4C{\gamma }_2{a}^2{S}^{(D-1)/D},(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

$\mathrm{\Phi }\approx S{a}^3⟨p⟩-{\gamma }_1{a}^2{S}_e.(34)$

If the invaded surface area S_ea² times the characteristic pore size a is proportional to the invaded volume Sa³, the ratio ϕ/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 ρ₁ and viscosity μ₁ = 0 invades another fluid with a higher density ρ₂ and viscosity μ₂ 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)d{p}_t.(35)$

We therefore have

$f-f_c={\int }_{p_c}^{p}N(p_t)d{p}_t,(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).(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(\mathrm{\Delta }\rho ga-\frac{q{\mu }_{2}a}{kA}\right)(x-x_0)/a=F(x-x_0)/a,(38)$

where x₀ corresponds to a the 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 Δρ = ρ₂ − ρ₁ the density difference between the fluids. The fluctuation number F

$F=N(p_c)\left(\mathrm{\Delta }\rho ga-\frac{q{\mu }_{2}a}{kA}\right),(39)$

is a dimensionless number [61] which characterizes the gravitational and viscous fields and the capillary pressure fluctuations. Note that the width of the capillary pressure threshold distribution is σ ≈ 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 η = (x − x₀) ∝ ξ, where ξ 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 ξ is given by [76].

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

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

$\eta /a\propto F^{\frac{-\nu }{1+\nu }}.(41)$

Therefore, when F > 0, the front will be stabilized with a characteristic length scale η [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 = Δρga²/γ = 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).

FIGURE 10

FIGURE 10. The dependence of the front width η/a on the generalized Bond number Fa/(γN (p_c)) = Bo − Ca. The blue dots correspond to data taken from [61] where the Bond number Bo = Δρga²/γ = 0.154 is kept fixed but the capillary number Ca = qμ₂a²/(γkA) is changed by varying the flow rate q.

Let us first assume that the gravitational effect is large enough such that it is the characteristic length scale η that sets the width of the front η < 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 η², 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,S_e\propto {\left(\frac{w}{\eta }\right)}^2{\left(\frac{\eta }{a}\right)}^{D_e},(42)$

and the work can be written as

$\begin{array}{l}\hfill W=-g({\rho }_{2}h+{\rho }_1(L-h))q{t}_w+\frac{q{\mu }_{2}h}{kA}S{a}^{2}b\\ \hfill +S{a}^{2}b⟨p⟩,\end{array}(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 Φ plus the surface energy E_s minus the change in gravitational potential energy ΔU_p (corresponding to the first term in Eq. 43)

$W=\mathrm{\Phi }+E_s-\mathrm{\Delta }{U}_p.(44)$

We then get the dissipation

$\begin{array}{l}\hfill \mathrm{\Phi }=\frac{q{\mu }_{2}h}{\kappa A}S{a}^{2}b+S{a}^{2}b⟨p⟩\\ \hfill -{\gamma }_2{a}^{2}2S-C{\gamma }_{2}ab2S^{(D-1)/D}-{\gamma }_{1}ab{S}_e^{D_e},\end{array}(45)$

which for large S can be approximated as

$\mathrm{\Phi }=\frac{q{\mu }_{2}h}{\kappa A}S{a}^{2}b+S{a}^{2}b⟨p⟩-{\gamma }_2{a}^{2}2S-{\gamma }_{1}ab{S}_e^{D_e}.(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 η > w, the relations for S and S_e will change to

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

but we will still have the same expressions for W and ϕ. Again if D_e = D, the ratios ϕ/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 η 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{L{w}^{d-1}}{{\eta }^d}\right){\left(\frac{\eta }{a}\right)}^D{a}^d}{L{w}^{d-1}}={\left(\frac{\eta }{a}\right)}^{D-d},(48)$

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

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

However, when the fluctuation number F is sufficient small, η > w, and the width of the model w will be the characteristic length scale in the problem. Then

$S_{nw}^F\propto {\left(\frac{w}{a}\right)}^{(D-d)}.(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 ≫ 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].

FIGURE 11

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μ₂a²/(γkA) = 1.2 ⋅ 10^–4 is kept fixed while the Bond number Bo = Δρga²/γ is changed in the experiments. Numerical and experimental data from [81].

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

$F=N(p_c)\left(\mathrm{\Delta }\rho ga+\frac{q{\mu }_{1}a}{k_{1}A}-\frac{q{\mu }_{2}a}{kA}\right),(51)$

where k₁ 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₁(F) will be a function of F. If we know this functional dependence we can solve the equation

$F=N(p_c)\left(\mathrm{\Delta }\rho ga+\frac{q{\mu }_{1}a}{k_1(F)A}-\frac{q{\mu }_{2}a}{kA}\right),(52)$

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

$\eta /a\propto F^{\frac{-\nu }{1+\nu }}.(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.

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.

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.

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.

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.

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

Google Scholar

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).

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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 Sourc (2020) 479:228708. doi:10.1016/j.jpowsour.2020.228708

CrossRef Full Text | Google Scholar

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

Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

61. Meheust 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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Hele-Shaw HS. The Flow of Water. Nature (1898) 58:34. doi:10.1038/058034a0

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar