Verification of a dynamic scaling for the pair correlation function during the slow drainage of a porous medium

TL;DR

This paper experimentally verifies a previously observed dynamic scaling law for the pair correlation function during slow drainage in porous media, connecting experimental data with theoretical models for different temporal regimes.

Contribution

It provides the first experimental validation of the dynamic scaling law for the pair correlation function in porous media drainage, linking theory and experiment.

Findings

Confirmed the scaling law $N(r,t) \,\propto\, r^{-1}f(r^{D}/t)$ experimentally.

Connected theoretical models to explain short-time and long-time behavior of $N(r,t)$.

Developed a new theoretical argument for the long-time regime exponent.

Abstract

We give experimental grounding for the remarkable observation made by Furuberg et al. in Ref. [furuberg1988] of an unusual dynamic scaling for the pair correlation function $N(r,t)$ during the slow drainage of a porous medium. The authors of that paper have used an invasion percolation algorithm to show numerically that the probability of invasion of a pore at a distance $r$ away and after a time $t$ from the invasion of another pore, scales as $N(r,t)\propto r^{-1}f\left(r^{D}/t\right)$, where $D$ is the fractal dimension of the invading cluster and the function $f(u)\propto u^{1.4}$, for $u \ll 1$ and $f(u)\propto u^{-0.6}$, for $u \gg 1$. Our experimental setup allows us to have full access to the spatiotemporal evolution of the invasion, which was used to directly verify this scaling. Additionally, we have connected two important theoretical contributions from the literature to explain the functional dependency of $N(r,t)$ and the scaling exponent for the short-time regime ($t \ll r^{D}$). A new theoretical argument was developed to explain the long-time regime exponent ($t \gg r^{D}$).

Equations

$$N(r,t)=r^{-1}f\left(\frac{r^D}{t}\right),$$

$$f(u)\propto \{\begin{array}{ll}{u}^{1.4}& \text{if u≪1}\\ u^{-0.6}& \text{if u≫1}\end{array}.$$

$$P_t(\Theta )\propto \frac{1}{t}{\left(\frac{\Theta }{t}\right)}^{-{\tau }_b^{all}}H\left(\Theta -t\right),$$

$$Q_{\Theta }(r)\propto {\Theta }^{-\left(1+\alpha \right)/D}{r}^{\alpha }H\left({\Theta }^{1/D}-r\right),$$

$$N(r,t)={\int }_0^{\infty }{P}_t(\Theta )Q_{\Theta }(r)d\Theta ,$$

$$N(r,t)r=f\left(\frac{r^D}{t}\right)=\{\begin{array}{ll}{\left(\frac{r^D}{t}\right)}^{\frac{1+\alpha }{D}}& \text{if trD≪1}\\ {\left(\frac{r^D}{t}\right)}^{1-{\tau }_b^{all}}& \text{if trD≫1}\end{array}.$$

$$1-{\tau }_b^{all}=-1+\left(\frac{D_e-1/\nu }{D}\right)⟹1-{\tau }_b^{all}=-0.68,$$

$$N(r,t)\propto \frac{1}{\left(t_1-t_{\ast }\right)}\frac{1}{t^{D_e/D}}\approx \frac{1}{t^{\xi }},$$

$$N(r,t)\propto \frac{1}{r}{\left(\frac{r^D}{t}\right)}^{\xi }.$$

Full text

Verification of a dynamic scaling for the pair correlation function during the slow drainage of a porous medium

M. Moura Corresponding author (marcel.moura@fys.uio.no) PoreLab, Department of Physics, University of Oslo, Norway

K. J. Måløy

PoreLab, Department of Physics, University of Oslo, Norway

E. G. Flekkøy

PoreLab, Department of Physics, University of Oslo, Norway

R. Toussaint

PoreLab, Department of Physics, University of Oslo, Norway

Université de Strasbourg, CNRS, IPGS UMR 7516, F-67000 Strasbourg, France

We give experimental grounding for the remarkable observation made by Furuberg et al. in Ref. [1] of an unusual dynamic scaling for the pair correlation function $N(r,t)$ during the slow drainage of a porous medium. The authors of that paper have used an invasion percolation algorithm to show numerically that the probability of invasion of a pore at a distance $r$ away and after a time $t$ from the invasion of another pore, scales as $N(r,t)\propto r^{-1}f\left(r^{D}/t\right)$, where $D$ is the fractal dimension of the invading cluster and the function $f(u)\propto u^{1.4}$, for $u\ll 1$ and $f(u)\propto u^{-0.6}$, for $u\gg 1$. Our experimental setup allows us to have full access to the spatiotemporal evolution of the invasion, which was used to directly verify this scaling. Additionally, we have connected two important theoretical contributions from the literature to explain the functional dependency of $N(r,t)$ and the scaling exponent for the short-time regime ($t\ll r^{D}$). A new theoretical argument was developed to explain the long-time regime exponent ($t\gg r^{D}$).

\saythanks

The slow drainage of a porous medium is characterized by a rich intermittent dynamics of invasion bursts, typically occurring at several time and length scales [2, 3, 4]. Similar intermittent activity is observed in a wide variety of physical, biological and social systems [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. The ubiquity of intermittent phenomena is an indication that its origin is not expected to depend on specific system-dependent details. It is generally associated to the competition between an adaptive external driving force and an internal random resistance against that force [3, 19]. Energy is slowly injected by the external force during stable periods which are abruptly interrupted by sudden dissipative events occurring at a much faster time scale [20, 21, 22]. In the case of two-phase flows in porous media, the external force could come from a syringe pump or an applied pressure difference across the sample, while the internal resistance is caused by the narrower or wider pore-throats that are invaded during the flow [2, 3, 4].

The balance between viscous, capillary and gravitational forces [23, 24, 25, 26] generates several interesting invasion patterns in porous media flows, ranging from compact invasion [27, 28, 29] to fractal structures [30, 31, 32, 33]. Although the pattern morphology has been well studied, very few studies have focused on its dynamical features [3, 4, 19, 34, 35], due in part to the difficulty of simultaneously analyzing detailed invasion data in both temporal and spatial domains. Single pore invasion events in rocks have only recently been imaged in real-time via modern X-ray microtomography techniques [34, 36], but extending those techniques to study pattern formations in rocks at the larger scales, while keeping single-pore resolution, remains a challenge.

The experimental constraints fueled the development of several numerical algorithms. Invasion percolation (IP) [37, 38], diffusion-limited aggregation (DLA) [39, 40, 30] and anti-DLA [40] have been employed in the simulation of, respectively, slow drainage (capillary fingering), fast drainage (viscous fingering) and stable imbibition (compact growth) [41]. In the standard IP model the invasion happens one pore at a time, always at the widest available pore-throat (with the smallest capillary pressure threshold). Using an IP model, Furuberg et al. [1] have addressed the following question: given that a reference pore located at position $\bm{r_{0}}$ is invaded at a time $t_{0}$, what is the probability that a second pore located at position $\bm{r_{1}}$ is invaded at some later time $t_{1}$? After the vanishing of transitional effects, the probability should only be a function of the relative quantities $r=\left|\bm{r_{1}}-\bm{r_{0}}\right|$ and $t=t_{1}-t_{0}$. The authors have defined a pair correlation function $N(r,t)$, such that $N(r,t)drdt$ gives the answer to the question, i.e., it is the probability of invasion of a pore located between distances $r$ and $r+dr$ and at a time between $t$ and $t+dt$ with respect to the invasion of the reference pore. By considering theoretical arguments related to the normalizaion of $N(r,t)$ and the connection between $N(r,t)$ and the pair connectedness function, the authors have suggested the form

[TABLE]

where the dynamic exponent $D$ corresponds to the fractal dimension of the invaded front [42, 43]. The IP simulations employed in that study numerically confirmed the validity of Eq. (1) and found out additionally that the function $f(u)$ presents the unusual scaling

[TABLE]

In the current work, we have demonstrated these results experimentally, almost 30 years after the original findings of Furuberg et al. [1].

Fig. 1 shows a diagram of the experimental setup. A quasi-2D porous network is formed by a monolayer of glass beads, with diameters $a$ in the range $1.0$\mathrm{m}\mathrm{m}$<a<1.2$\mathrm{m}\mathrm{m}, randomly placed in the gap of a modified Hele-Shaw cell. The beads are kept in place by a pressurized cushion on the bottom plate of the cell. The porous medium is initially saturated with a viscous liquid (wetting phase), composed of glycerol ($80\%$ in weight) and water ($20\%$ in weight), having kinematic viscosity $\nu=4.25\cdot 10^{-5}\,$\mathrm{m}^{2}\mathrm{/}\mathrm{s}, density $\rho=1.205\,$\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3} and surface tension $\gamma=0.064\,$\mathrm{N}\mathrm{/}\mathrm{m} (see Ref. [44] for additional details). The outlet channel is connected to a syringe pump, from which liquid can be slowly withdrawn at a constant rate of $q=0.0050\,$\mathrm{m}\mathrm{l}\mathrm{/}\mathrm{m}\mathrm{i}\mathrm{n}$\approx 0.14\,$\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{/}\mathrm{s}, thus assuring that the dynamics happens in the capillary regime, in which capillary forces dominate over viscous ones [[33](#bib.bib33)]. Air (non-wetting phase) enters from a width-spanning inlet channel, which is kept open to the atmosphere. Once the capillary pressure (difference in pressure between the non-wetting and wetting phases) is large enough to overcome the threshold associated with a given pore-throat, the invasion of one or more pores happens, and viscous pressure drops are triggered within the medium, thus dissipating energy [[3](#bib.bib3), [45](#bib.bib45)]. Pictures are taken every $34\,$\mathrm{s} which leads to an average number of $K\approx 4.7$ invaded pores per image (average pore invasion time of $t_{p}=7.2\,$\mathrm{s}). The capillary number is, $Ca=\rho\nu q/\Sigma\gamma=6.1\cdot 10^{-7}$, where $\Sigma=1.1\cdot 10^{-4}\,$\mathrm{m}^{2} is the cross section of the model. The choice of a slow constant withdrawal rate reflects the fact that in the IP algorithm employed in Ref. [1], viscous pressure drops are neglected and exactly 1 pore (numerical site) is invaded per time step.

Fig. 2 shows the spatiotemporal map of the invasion. The flow is from left to right but statistically similar results should be found if the flow were in the opposite direction since the model is prepared randomly. The experiment stops at breakthrough, i.e., when the air phase first percolates (equivalent to the infinite cluster in IP). Only the central $90\%$ of the model’s length (inlet–outlet direction) is considered, in order to avoid boundary effects which have been observed close to the inlet and outlet [44]. The information content in this map is similar to that in the IP simulations, namely the position and invasion times of all pores, thus allowing the computation of the pair correlation function $N(r,t)$. In order to precisely locate the pores, we performed a Delaunay triangulation [46] of the points marking the centers of all glass beads and then identified the centroid of each Delaunay triangle as a pore center. In the IP simulations [1], the time $t$ is naturally measured as the number of invaded pores, whereas in our experiments it is given by the image number $t_{I}$. The conversion between $t_{I}$ and $t$ is $t_{I}=t/K$, where $K$ is the average number of pores invaded per image. This linear conversion works best for large times (several pore invasions) since the intermittency of the Haines jumps [2] is relevant for short times.

For a given reference pore, one can produce a histogram of the euclidean distances $r$ of all pores invaded between times $t_{I}$ and $t_{I}+\Delta t_{I}$ after its own invasion. By treating each invaded pore in the system as a reference for the production of one such histogram and adding them up, we obtain the function $N(r,t)$ (apart from a normalization factor). We have chosen $\Delta t_{I}=2$ images to guarantee that a small number of pore invasions (given on average by $\Delta t_{I}K\approx 9.4$) would be present on each histogram. Different values of $\Delta t_{I}$ in the range $1\leq\Delta t_{I}\leq 5$ were tested and the results presented did not change significantly, since the experiment duration ($\approx 33\,$\mathrm{h}) is much longer than the time between images ($\approx 34\,$\mathrm{s}). It is worth mentioning that, due to the fast nature of the Haines jumps [2], if one wants to capture the exact invasion time of a pore, one would need to rely on high-speed imaging, which would then dramatically constrain the total duration of the experiment. Such level of accuracy was not needed in this study.

The measured function $N(r,t)$ is shown in Fig. 3 for $9$ fixed values of time $t/K=t_{I}$ (shown in the legend). The distance $r$ is measured in pore length units (euclidean distance divided by a characteristic interpore distance $r_{p}=0.92\,$\mathrm{m}\mathrm{m}$$). $N(r,t)$ presents a peak indicating a maximum probability of invasion at a certain relative distance $r_{t}$ which increases with time (since the air-liquid interface had more time to move, pores farther away can be invaded). The validity of Eq. (1) is proved by using this equation to collapse the data from Fig. 3. The result, shown in Fig. 4, indicates that the product $N(r,t)\>r$ is indeed a function only of the reduced variable $u=r^{D}/t$ and not of $r$ and $t$ separately (the fractal dimension of the invasion cluster was measured to be $D=1.75$). This confirms the validity of the dynamic scaling in Eq. (1). The fact that the product $N(r,t)\>r$ is peaked around $r^{D}/t=1$ (with $r$ and $t$ measured in terms of pore length and average pore invasion time units) indicates that the most probable place for invasion occurs at a distance $r=t^{1/D}$, as noted in [1]. If one chooses another set of measuring units, the maximum changes to $r_{p}^{D}/t_{p}$, where $r_{p}$ and $t_{p}$ are the typical interpore distance and pore invasion time in the new units. The scaling behavior of $N(r,t)$ could not be easily inferred from the separate curves in Fig. 3, due to the limited statistics, but becomes more visible after the data collapse in Fig. 4. Both scaling regimes, $u^{1.4}$ for $u\ll 1$ and $u^{-0.6}$ for $u\gg 1$, are well reproduced by the experiments (guide-to-eye thick solid lines). Indeed, even the deviation from the scaling observed as the dropping curves for large times and $u\gg 1$ (due to finite-size effects) are also in agreement with the simulations in [1].

Next we analyze the origin of this scaling behavior. In the work of Roux and Guyon [47], the distribution of temporal avalanches in the time series of threshold values from an IP model was used to propose an analytical prediction for the unusual scaling presented in Eqs. (1) and (2). Later on, Maslov [48] pointed out some inconsistencies in the assumptions used in Ref. [47], which interfered with their estimation for the exponents in Eq. (2). Nevertheless, the reasoning in Ref. [47] to justify the functional form behind Eqs. (1) and (2) still applies. By taking into account the results of Maslov [48], into the analysis performed by Roux and Guyon [47], we obtain a consistent value for the short-time exponent in Eq. (2) (for $u\gg 1$ or $t\ll r^{D}$). The analysis for the long-time exponent (for $u\ll 1$ or $t\gg r^{D}$) will require some additional considerations, as shown later.

Following Ref. [47], for a given time interval $t$, the distribution of backwards temporal avalanches $\Theta$ (as defined in [47] and also in [48]) that are larger than $t$ is assumed to follow the power law

[TABLE]

where $H(x)=1$ for $x>0$ and $H(x)=0$ otherwise. The dependency on $t$ in the prefactor follows from the normalization in the interval $t<\Theta<\infty$ [47]. We remind that in the IP model, the time corresponds to the size (mass) of the invaded cluster, measured in number of pores. For a cluster of size $\Theta$, the distribution $Q_{\Theta}(r)$ of distances $r$ between pores in that cluster is also assumed to follow a power law,

[TABLE]

for $r$ in the interval $0<r<\Theta^{1/D}$, where again the dependency on $\Theta$ in the prefactor is obtained by the normalization in the interval $0<r<\Theta^{1/D}$ [47]. The pair correlation function $N(r,t)$ is then given by

[TABLE]

which, using Eqs. (3) and (4) leads to

[TABLE]

The insightful argumentation provided by Roux and Guyon [47] naturally generates the correct functional form of Eqs. (1) and (2). Nevertheless their analysis for the numerical values of the exponents presents some inconsistent assumptions, as noted by Maslov [48]. The exponent $\tau_{b}^{all}$ characterizing the distribution of all backwards avalanches was shown [48] to be given by $\tau_{b}^{all}=3-\tau$, where $\tau=1+\left(D_{e}-1/\nu\right)/D$ is the cluster size distribution exponent (first derived in [49], and verified experimentally in [50]), $D$ and $D_{e}$ are respectively the fractal dimensions of the growing cluster and its boundary and $\nu=4/3$ is the exponent characterizing the divergence of the correlation length [43, 38]. Using the values $D=1.82$ and $D_{e}=4/3$ (external perimeter growth) [43, 38], we find that the short-time exponent is

[TABLE]

which is consistent with our measurements and very close to the value $-0.6$ reported in Ref. [1], see Fig. 4.

The computation of the long-time exponent brings additional challenges due to the difficulty in estimating the exponent $\alpha$ in Eq. (4) [47]. With that in mind, we take an alternative approach here. Consider the situation in which a pore-throat that gives access to a pore at position $\bm{r_{1}}$ has been reached by the liquid air interface at a time $t_{*}$. From that moment on, that pore-throat and the pore-body at $\bm{r_{1}}$ are available to the invasion. Let us focus here on the case in which this particular pore-throat has a relatively high value of capillary pressure threshold, such that its invasion will typically take a long time to occur. This type of event lies in the long-time regime, for which $t\gg r^{D}$, and we are interested in addressing their relative probability of occurrence.

For the invasion of the pore at $\bm{r_{1}}$ to happen exactly at time $t=t_{1}$, a set of two well defined conditions must be verified at that time: 1) the capillary pressure must reach a historically high value since time $t=t_{*}$ and 2) the pore-throat at $\bm{r_{1}}$ must have the lowest value of capillary pressure threshold (i.e., be the widest pore-throat) among those that belong to the liquid air interface. Let us initially address condition 1). Consider the capillary pressure signal as the discrete time series formed by the sequence of the capillary pressure thresholds $p\left(t\right)$ associated with the invasion of successive pore-throats. If the system is in a statistical steady state (say, the flow has been going for a long time and happens in a long rectangular cell of large but finite width), $p\left(t\right)$ fluctuates around some well defined average and the historical maximum is equally likely to occur anywhere in the interval $t_{*}<t\leq t_{1}$. The probability that it occurs at the extreme point $t=t_{1}$ is then proportional to $1/\left(t_{1}-t_{*}\right)$. Once the capillary pressure reaches the historical maximum at $t=t_{1}$, the invasion at $\bm{r_{1}}$ can happen. Next, we consider condition 2) and calculate the probability that the particular pore-throat in question is the one with the lowest threshold. Since the capillary pressure at time $t=t_{1}$ has reached a historical maximum, it is the first time that the pore-throats on the interface have been tested against such a high value of capillary pressure, and at this pressure they all have the same invasion probability. The number of sites $N_{I}$ that belong to the interface of the cluster that has grown since time $t=t_{0}$ (when the reference pore at $\bm{r_{0}}$ was invaded) scales as $N_{I}\propto t^{D_{e}/D}$. The probability that the particular pore-throat at $\bm{r_{1}}$ is the widest is simply given by $1/N_{I}$. (We have made the assumption that the invasion of the pore at $\bm{r_{1}}$ does not depend on invasion events that happened earlier than the invasion at $\bm{r_{0}}$. In the limit $t=t_{1}-t_{0}\to\infty$, this approximation becomes exact). By considering the probability of simultaneously satisfying conditions 1) and 2), we have

[TABLE]

where $\xi=1+D_{e}/D$ and, in the limit of large times, $t_{1}-t_{*}\approx t_{1}-t_{0}=t$. The $r$ dependence in the previous equation is obtained by considering again Eq. (1). Since the product $N(r,t)r$ is a function only of the reduced variable $r^{D}/t$, we have that in the long-time regime,

[TABLE]

Using the literature values $D=1.82$ and $D_{e}=4/3$ [43, 38], we find $\xi=1.73$. This value is higher than the exponent $\xi=1.4$ reported by Furuberg et al. [1], but is consistent with the experimental data, particularly in the extreme left region of Fig. 4, where $t\gg r^{D}$ and the approximations of the model should hold best.

After the aged (hard) site at $\bm{r_{1}}$ is invaded, we could in principle have the invasion of easier pores in its vicinity that would also contribute to the counting in the long-time regime of $N(r,t)$. Although the argument made here counts explicitly only the aged (hard) sites, the invasion of easier ones in the long-time regime is conditioned to the prior invasion of an aged site and therefore should not change the temporal scaling of Eq. (8).

In this work we have given experimental validation to the unusual dynamic scaling for the pair correlation function $N(r,t)$ during the slow drainage of a porous medium, first observed by Furuberg et al. [1] nearly 30 years ago. Although this important result has been reproduced in other numerical works [51], to the best of our knowledge, the experimental verification presented here is new. By linking two important contributions to the literature, the works of Roux and Guyon [47] and Maslov [48], we found out that both approaches lead to the same predictions for the short-time exponent of $N(r,t)$, which agrees well with our measurements. A new theoretical explanation for the long-time exponent has also been provided, with good agreement with the experimental data. Possible extensions of the work include other flow regimes and geometries (e.g. radial injection). While $N(r,t)$ probably varies for faster flows, it possibly remains unchanged for slow drainage in a 2D radial geometry, once transient effects from the point inlet are dissipated.

We gratefully acknowledge the support from the University of Oslo, University of Strasbourg, the Research Council of Norway through its Centre of Excellence funding scheme with project number 262644, the CNRS-INSU ALEAS program and the EU Marie Curie ITN FLOWTRANS network.