ORIGINAL RESEARCH article
Burst Dynamics, Upscaling and Dissipation of Slow Drainage in Porous Media
- 1PoreLab, Department of Physics, The Njord Centre, University of Oslo, Oslo, Norway
- 2PoreLab, Department of Geoscience and Petroleum, Norwegian University of Science and Technology, Trondheim, Norway
- 3PoreLab, Department of Physics, Norwegian University of Science and Technology, Trondheim, Norway
- 4PoreLab, Department of Chemistry, Norwegian University of Science and Technology, Trondheim, Norway
- 5CNRS, 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.
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  and solar cell technology , fiber-reinforced composite materials , textile fabric characterization , 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 . 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  a topic of relevance for the development of modern battery technology .
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  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 . 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 , 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  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 Es 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, Es, 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 Es/W and Φ/W are both independent of system size.
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 = g0 sin(θ) along the x direction of the model, where g0 = 9.82 m/s2. 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 , 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  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 m2/s, ρ = 1.205 g/cm3 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 r0 to one with a radius of curvature r = r0—dr. This results in a volume change dvi of the nonwetting fluid, and a corresponding volume change—dvi 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 r0 to the front defined by the radius of curvature r = r0—dr. The longest of these distances we denote as dxi. The distance dxi to lowest order in dr is given by
depends on the wetting angle ψ (See Figure 2) and the capillary pressure p0 given by the Young-Laplace equation p = γ(1/r + 1/r1)  for r = r0. Here r is the in-plane, and r1 the out of plane radius of curvature assumed to be constant in the quasi-two-dimensional experiments.
FIGURE 2. The figure illustrates a nonwetting fluid invading a wetting fluid at a single pore-throat at different radius of curvature r0 and r = r0 − dr. The distance moved between the two front positions is dxi and ψ is the wetting angle.
We can then calculate the corresponding increase in capillary pressure
The volume change dvi of the invading fluid in the pore can be written as
where the capacitive volume
gives the fluid volumetric change per unit capillary pressure in pore i .
Let Ui(x) be the elastic energy (surface energy) of the interface in pore i, and dxi 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 xi would then be
We have assumed p0 = 0 without loss of generality since
where κi is the capacitive volume of pore i so that κ is the average capacitive volume over the n interface pores. The capillary pressure p2 at the time t2 and the capillary pressure p1 at the time t1 are related by
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) . The capillary pressure will then build up from p1 at t1 to p2 at t2 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 Ut (t1) to Ut (t2)
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 pt. 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 t2 at a capillary pressure p2 = pt and that, after the burst, the system reaches another equilibrium state at capillary pressure p3. 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
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 = p2 − p3 is the pressure drop in the burst and (p2 + p3)/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γ2a2s and one from the interface contour aγ1bse. 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 ase is the length of the invasion contour of the burst (including trapped clusters). The factor 2 in the term 2γ2a2s 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, γ1 and γ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), De of the invasion contours, and Df of the invasion front are shown in Figure 4B), where a box counting technique was employed in the measurement . We have obtained the values D = 1.87 ± 0.10, De = 1.84 ± 0.10 and Df = 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. (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 De of the invasion contours and the fractal dimension Df 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 γ1 = γ2 = γ which is the surface tension between the two fluids. However when no film is left behind the invasion front, γ2 = γ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 γ2 = γ 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 ϵ = Ans/At, where Ans is the interface’s surface area between the nonwetting phase and the solid phase, and At 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 = (1 − ϵ)γ + ϵ(γns − γws) for the total internal interfacial contours. Using again the Young equation (10) we can rewrite this expression as γ1 = (1 − ϵ)γ + ϵγ cos(ψ).
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 Es, such that
It is important to note that W, Φ, and Es are average quantities that do not account for fluctuations, and that ϕ and Es are averaged over a large number of bursts.
From  we know that the distribution function of the burst sizes G(s) is directly linked to the distribution of the capillary pressure drops
where the cutoff size of the bursts s∗ is .
where N (pc) is the value of the normalized capillary pressure threshold distribution taken at the critical percolation threshold pressure pc, and ν = 4/3 is the exponent describing the divergence of the correlation length in percolation . Since the capillary pressure threshold distribution N (pt) is normalized, the width of the distribution σ ≈ 1/N (pc). Therefore the cutoff size of the bursts will increase with σ according to Eq. 20 as s∗ ∝ σ0.78. Martys et al.  derived the analytical form
where D and Df 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 Df = 4/3 [23, 27, 58], we obtain τ = 1.32, consistent with our experimental data shown in 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. .
FIGURE 7. The dependence of the crossover function H (Δp/Δp*) of the pressure jump distribution in a modified invasion percolation simulation . 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 .
We will in the following use
where M is the total number of bursts during the time considered. Here we have assumed that the average pressure in each burst (p2 + p3)/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. 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.
where I is the integral
Then by using Eq. 23 we obtain
Assume that we consider the time tw 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 a2S 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
where I is the integral in (26). We then get the following simple expression for the work
This expression is the well-known relation of the work expressed as a pressure times a volume change ⟨p⟩ΔV, where
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 Se seen in red in Figure 4. For large S the third term can be neglected and
The length of the front n, seen as the green line in Figure 4, is given by
where Df is the fractal dimension of the front (not included trapping) and De is the fractal dimension of contours of the trapped clusters. From Figure 4 we see that the measured fractal dimensions for De and D seem very close, effectively indistinguishable from one another considered the error bars. This result implies that S ∝ Se. In the experiment shown in Figure 4 we have measured the proportionality relation to be Se = 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 Est = γ2a22S (green) and the lateral interface Esl = γ1abSe (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, Esl/W and Est/W are independent of w/a for large systems.
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 = (1 − ϵ)γ + ϵγ cos(ψ) = 0.034N/m, ⟨p⟩ = 210Pa and D = 1.87, see Figure 4. As previously stated, our assumptions for γ1 and γ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 w2, Eq. 29 for the work would remain the same, while Eq. 30 for the dissipation would be modified to
where C is a constant. When the system gets large the last term can be neglected compared to the two first terms such that
If the invaded surface area Sea2 times the characteristic pore size a is proportional to the invaded volume Sa3, the ratio ϕ/W and Esl/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 ρ1 and viscosity μ1 = 0 invades another fluid with a higher density ρ2 and viscosity μ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].
We therefore have
where the critical occupation probability is fc and the critical percolation pressure is pc. As before, N (pt) is the normalized capillary pressure threshold distribution. Now, by Taylor expanding N (pt) around pc to the lowest order in p − pc in Eq. 36 we get
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
where x0 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 Δρ = ρ2 − ρ1 the density difference between the fluids. The fluctuation number F
is a dimensionless number  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 (pc), 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 − x0) ∝ ξ, 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 .
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  where the authors analyzed the scaling of the invasion front in drainage by keeping the Bond number Bo = Δρga2/γ = 0.154 fixed but changing the capillary number Ca by varying the flow rate q (see Figure 10 in Ref.  and experiments 1, 2, 3, 4, 5, 7, and 9 in Table I in the same reference).
FIGURE 10. The dependence of the front width η/a on the generalized Bond number Fa/(γN (pc)) = Bo − Ca. The blue dots correspond to data taken from  where the Bond number Bo = Δρga2/γ = 0.154 is kept fixed but the capillary number Ca = qμ2a2/(γ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 tw corresponding to the time the front needs to move a length scale w. Dividing space into cells of size η2, below which the fluid distribution is fractal, the number of pores S and Se are given by
and the work can be written as
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 Es minus the change in gravitational potential energy ΔUp (corresponding to the first term in Eq. 43)
We then get the dissipation
which for large S can be approximated as
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 Se will change to
but we will still have the same expressions for W and ϕ. Again if De = D, the ratios ϕ/W, Esl/W and Est/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. . 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
where d is the spatial dimension (2 or 3). Hence using Eq. 41
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
Figure 11 shows a two-dimensional invasion percolation simulation with a gravitational field together with drainage experiments performed by Ayaz et al. . 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 .
FIGURE 11. The final saturation as a function of the generalized Bond number Fa/(γN (pc)) = 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μ2a2/(γkA) = 1.2 ⋅ 10–4 is kept fixed while the Bond number Bo = Δρga2/γ is changed in the experiments. Numerical and experimental data from .
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
where k1 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 k1(F) will be a function of F. If we know this functional dependence we can solve the equation
with respect to F and thereby find the characteristic length scale η from
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 η.
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 fc. The work W, dissipation Φ, and surface energy Es 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 . 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.
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.
This work was partly supported by the Research Council of Norway through its Centers of Excellence funding scheme, project number 262 644.
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.
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.
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
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
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).
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
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
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
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
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
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
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
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
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 Sourc (2020) 479:228708. doi:10.1016/j.jpowsour.2020.228708
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
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
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
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
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
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
Keywords: porous media, burst dynamics, intermittency, drainage, multiphase flow, critical phenomena, percolation
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
Received: 15 October 2021; Accepted: 15 November 2021;
Published: 10 December 2021.
Edited by:Fernando A. Oliveira, University of Brasilia, Brazil
Reviewed by:Boqi Xiao, Wuhan Institute of Technology, China
Hans Herrmann, UMR7636 Laboratoire de physique et mécanique des milieux hétérogenes (PMMH), France
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.
*Correspondence: Knut Jørgen Måløy, firstname.lastname@example.org