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

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

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 [1][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 systemdependent details. It is generally associated with 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 dynamic 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 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 [30,39,40], and antidiffusion-limited aggregation [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 r 0 is invaded at a time t 0 , what is the probability that a second pore located at position 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 ¼ jr 1 − r 0 j 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 normalization of Nðr; tÞ and the connection between Nðr; tÞ and the pair connectedness function, the authors have suggested the form 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 In the current Letter, we have demonstrated these results experimentally, almost 30 years after the original findings of Furuberg et al. [1]. Figure 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 < a < 1.2 mm, 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% by weight) and water (20% by weight), having kinematic viscosity ν ¼ 4.25 × 10 −5 m 2 =s, density ρ ¼ 1.205 g=cm 3 , and surface tension γ ¼ 0.064 N=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 ml= min ≈0.14 pore=s, thus assuring that the dynamics happens in the capillary regime, in which capillary forces dominate over viscous ones [33]. Air (nonwetting phase) enters from a width-spanning inlet channel, which is kept open to the atmosphere. Once the capillary pressure (difference in pressure between the nonwetting 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,45]. Pictures are taken every 34 s, which leads to an average number of K ≈ 4.7 invaded pores per image (average pore invasion time of 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 one pore (numerical site) is invaded per time step. Figure 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 þ Δ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 Δt I ¼ 2 images to guarantee that a small number of pore invasions (given on average by Δt I K ≈ 9.4) would be present on each histogram. Different values of Δt I in the range 1 ≤ Δt I ≤ 5 were tested and the results presented did not change significantly, since the experiment duration (≈33 h) is much longer than the time between images (≈34 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 a level of accuracy was not needed in this study.
The measured function Nðr; tÞ is shown in Fig. 3 for nine 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 mm). 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 Ref. [1]. If one chooses another set of measuring units, the maximum changes to r D p =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 ≪ 1 and u −0.6 for u ≫ 1, are well reproduced by the experiments (guide-to-eye thick solid lines). Indeed, even the deviations from the scaling observed as the dropping curves for large times and u ≫ 1 (due to finite-size effects) are also in agreement with the simulations in Ref. [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] in the analysis performed by Roux and Guyon [47], we obtain a consistent value for the short-time exponent in Eq. (2) (for u ≫ 1 or t ≪ r D ). The analysis for the long-time exponent (for u ≪ 1 or t ≫ 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 Θ (as defined in Ref. [47] and also in Ref. [48]) that are larger than t is assumed to follow the power law 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 < Θ < ∞ [47]. Recall 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 Θ, the distribution Q Θ ðrÞ of distances r between pores in that cluster is also assumed to follow a power law, for r in the interval 0 < r < Θ 1=D , where again the dependency on Θ in the prefactor is obtained by the normalization in the interval 0 < r < Θ 1=D [47]. The pair correlation function Nðr; tÞ is then given by which, using Eqs. (3) and (4), leads to 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 τ all b characterizing the distribution of all backwards avalanches was shown [48] to be given by is the cluster size distribution exponent (first derived in Ref. [49], and verified experimentally in Ref. [50]), D and D e are, respectively, the fractal dimensions of the growing cluster and its boundary, and ν ¼ 4=3 is the exponent characterizing the divergence of the correlation length [38,43]. Using the values D ¼ 1.82 and D e ¼ 4=3 (external perimeter growth) [38,43], we find that the short-time exponent is 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 α 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 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 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 ≫ r D , and we are interested in addressing its relative probability of occurrence.
For the invasion of the pore at 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 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ðtÞ 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ðtÞ fluctuates around some well defined average and the historical maximum is equally likely to occur anywhere in the interval t Ã < t ≤ t 1 . The probability that it occurs at the extreme point t ¼ t 1 is then proportional to 1=ðt 1 − t Ã Þ. Once the capillary pressure reaches the historical maximum at t ¼ t 1 , the invasion at 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 r 0 was invaded) scales as N I ∝ t D e =D . The probability that the particular pore throat at r 1 is the widest is simply given by 1=N I . (We have made the assumption that the invasion of the pore at r 1 does not depend on invasion events that happened earlier than the invasion at r 0 . In the limit t ¼ t 1 − t 0 → ∞, this approximation becomes exact). By considering the probability of simultaneously satisfying conditions (1) and (2), we have where ξ ¼ 1 þ D e =D and, in the limit of large times,  Fig. 4, where t ≫ r D and the approximations of the model should hold best. After the aged (hard) site at 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 Letter 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 No. 262644, the CNRS-INSU ALEAS program, and the EU Marie Curie ITN FLOWTRANS network.