Critical behavior in porous media flow

The intermittent burst dynamics during the slow drainage of a porous medium is studied experimentally. We have shown that this system satisfies a set of conditions known to be true for critical systems, such as intermittent activity with bursts extending over several time and length scales, self-similar macroscopic fractal structure and $1/f^\alpha$ power spectrum. Additionally, we have verified a theoretically predicted scaling for the burst size distribution, previously assessed via numerical simulations. The observation of $1/f^\alpha$ power spectra is new for porous media flows and, for specific boundary conditions, we notice the occurrence of a transition from $1/f$ to $1/f^2$ scaling. An analytically integrable mathematical framework was employed to explain this behavior.


Introduction
The topic of fluid motion inside a porous network has deservedly been subjected to a considerable number of studies over the past decades. Scientists have studied the morphology and dynamics of the flow [1][2][3][4][5][6][7][8][9][10][11][12] and proposed a set of numerical schemes able to reproduce the observed macroscopic patterns [13][14][15][16][17] and relevant pore-scale mechanisms [18][19][20][21][22][23][24][25][26]. The topic is also of central importance for the study of groundwater flows and soil contaminants treatment [27,28] and has direct applications in the energy sector, for example, in hydrocarbon recovery methods [29]. One particularly interesting aspect of multiphase flow in porous media is its intermittent dynamics [3,4,18], with long intervals of stagnation followed by short intervals of strong activity. This kind of general behavior [30][31][32] appears in many physical, biological and economical systems, such as the stick-slip motion of a block on an inclined plane [33], the propagation of a fracture front in a disordered material [34][35][36], the number of mutations in models of biological evolution [37], acoustic emissions from fracturing [38,39], variations in stock markets [40], and the rate of energy transfer between scales in fully developed turbulence [41,42]. Intermittent phenomena arise irrespectively of the (certainly different) specific details of each system. In the particular case of porous media flows, this is caused by the interplay between an external load (for example, an imposed pressure difference across the system) and the internal random resistance due to the broader or narrower pore-throats. * Corresponding author (marcel.moura@fys.uio.no) In the present work we show experimental results on the burst dynamics during drainage in artificial porous media and investigate the question of how the pressure fluctuations (due to the burst activity) can encode useful information about the system. The flows studied are slow enough to be in the capillary regime, in which capillary forces are typically much stronger than viscous ones [3,43]. We have employed synthetic quasi-2D systems driven by a controlled imposed pressure (CIP) boundary condition. This boundary condition differs from the controlled withdrawal rate (CWR), more commonly used [3,9]. The dynamics is characterized both via direct imaging of the flow and by local pressure measurements. We present results related to the statistics of bursts, their morphology and orientation within the medium, and the power spectral density (PSD) associated with the fluctuations in the measured pressure signal. In particular, we show that for systems driven by the CIP boundary condition, the PSD presents a 1/f scaling regime. The presence of 1/f α power spectra is a widespread feature occurring in a myriad of contexts [44][45][46], commonly signaling the collective dynamics of critical systems. Some examples are the early measurements of flicker noise in vacuum tubes [47], fluctuations in neuronal activity in the brain [48], quantum dots fluorescence [49], loudness in music and speech [50,51] and fluctuations in the interplanetary magnetic field [52]. Although 1/f α power spectra have also been observed in some fluid systems, such as simulations and experiments on hydrodynamic and magnetohydrodynamic turbulence [53,54] and quasi-2D turbulence in electromagnetically forced flows [55], to the best of our knowledge the results reported here provide the first experimental observations of 1/f α power spectra in porous media flows. Fig. 1 shows a schematic representation of the setup employed (additional details in Ref. [56]). The quasi-2D porous network is formed by a modified Hele-Shaw cell filled with a monolayer of glass beads having diameters a in the range 1.0mm < a < 1.2mm. The beads are kept in place by a pressurized cushion placed on the bottom plate of the cell. A spongeous filter with pores much smaller than those in the medium is placed between the porous network and the outlet of the model. This filter allows the dynamics to continue inside the medium even after breakthrough [56]. Pressure measurements are taken at the outlet with an electronic pressure sensor (Honeywell 26PCAFG6G) that records the difference between the air pressure (non-wetting phase) and the liquid pressure (wetting phase) at the outlet, i.e., p m = p nw − p out w . Since the inlet is open to the atmosphere, p nw = p 0 in all experiments, where p 0 is the atmospheric pressure. The porous matrix was initially filled with a mixture of glycerol (80% in weight) and water (20% in weight) having kinematic viscosity ν = 4.25 10 −5 m 2 /s, density ρ = 1.205 g/cm 3 and surface tension γ = 0.064 N.m −1 . We have performed experiments on 4 different porous media with dimensions: (1) 27.3cm x 11.0cm, (2) 14.0cm x 11.5cm, (3) 32.8cm x 14.6cm and (4) 32.0cm x 4.5cm, where the first number corresponds to the length (inletoutlet direction) and the second to the width. The outlet of the model is connected to an external reservoir. The height difference h between the surface of the liquid in this reservoir and the model is used to control the imposed pressure via an adaptive feedback mechanism (CIP boundary condition). This mechanism guarantees that the pressure is only increased  The vast blue areas in the top image contain many smaller bursts (detail).

Methodology
when the system is in a quasi-equilibrium situation (see details in [56]). By slowly increasing the imposed pressure (via small steps in the height of the reservoir dh = 10µm =⇒ dp = ρgdh = 0.12 P a, where g is the acceleration of gravity), new pore-throats may become available to invasion. The value of dh was chosen to satisfy the accuracy condition that the height would typically have to be increased several times before new pores are invaded. As long as this condition is satisfied, the results obtained should be independent of the particular value of dh.

Burst size distribution
We begin by analyzing the size distribution of invasion bursts in a CIP experiment. A burst is understood as any connected set of pores invaded in the interval Θ = t 2 − t 1 between two consecutive time instants, t 1 and t 2 , at which the imposed pressure was increased (i.e., the imposed pressure is constant during the interval Θ, being changed only at its extremes t 1 and t 2 ). Fig. 2 shows the individual bursts for experiment CIP-1 (the number identifies the model), colored according to their area (top) and randomly (bottom), the latter being done to aid the visualization of separate bursts. Only bursts having their centroids in the central 90% of the length are considered, to avoid possible boundary effects [56]. A great deal of information can be obtained from this image. Initially, one can observe the homogeneity and isotropicality of the dynamics: the bursts don't seem to follow a well defined size gradient (the top image does not seem to transition from blue to red following a specific di- The line shows the scaling N (n) ∝ n −τ , with τ = 1.37 ± 0.08, which is consistent with the theoretical value τ = 1.30 ± 0.05 predicted by numerical simulations and percolation theory [21,58]. The data has been shifted vertically to aid visualization. rection), nor have they a clear preferred orientation (they are not particularly elongated in any direction). It is hard, if not impossible, to say from this image in which direction the invasion takes place (it is from left to right). A reflection (vertical or horizontal) or a 180 • rotation would also not be clearly identified. The box counting fractal dimension [57,58] of the invading cluster was measured to be D = 1.76 ± 0.05. Fig. 3 shows the burst size distribution N (n) for 3 separate experiments (the number of pores n being measured by normalizing the burst area by a typical pore area ≈ 0.3mm 2 ). The system exhibits the scaling N (n) ∝ n −τ , with τ = 1.37 ± 0.08, over at least two decades. The burst dynamics is therefore spatially self-similar, a feature commonly associated with systems close to a critical transition [46,57]. The exponent τ has been calculated via maximum likelihood estimation (MLE) [59] using the data from Fig. 2 for burst sizes in the interval 1 pore < n < 150 pores. MLE was used in order to avoid possible biases from data binning (MLE is a binning-free method), see also [60]. The scaling is shown in Fig. 3 on top of the logarithmically binned histogram of the data for the sake of visualization. Experiment CIP-4 was left out of the analysis because boundary effects rendered the results unreliable (model 4 is too narrow). The measured exponent is consistent with the value τ = 1.30 ± 0.05 predicted by numerical simulations and percolation theory [21,58]. Martys et al. [21] derived the analytical form where D and D e are respectively the fractal dimensions of the growing cluster and its external perimeter and ν ′ = 4/3 is the exponent characterizing the divergence of the correlation length [57,58]. Using the values D = 1.76 and D e = 4/3 [14], we obtain τ = 1.33, very close to the measured value τ = 1.37±0.08 shown in Fig. 3. Our measurements provide a direct experi- Crandall et al. [61] performed measurements in a CWR system finding the exponent τ = 1.53, which is compared to the theoretical prediction of τ = 1.527 from Roux and Guyon [62]. Nevertheless, Maslov [63] pointed out an inconsistency in this theoretical prediction, the correct expression being given in Eq. (1). Modified invasion percolation simulations and pressure measurements [3,4] have shown that, in a CWR, system very large bursts are split into smaller ones. A burst size distribution was observed, with exponent τ = 1.3 ± 0.05 for the simulations and τ = 1.45 ± 0.10 for the experiments (consistent with Eq. (1)), followed by an exponential cutoff [3,4]. In the CIP case large bursts can occur because the displaced liquid can freely flow out of the model but in the CWR case this is not possible since the available volume for the displaced liquid is bounded by the outlet syringe volume.

Burst time distribution
Let us now focus on the distribution G(Θ) of time intervals Θ between two successive increments in the imposed pressure during which invasion bursts have occurred. Fig. 4 shows this distribution, produced for all bursts with Θ > 120s, a cutoff related to the minimum time difference for proceeding the image analysis used in the feedback mechanism [56]. It scales as G(Θ) ∝ Θ −γ with γ = 2.04 ± 0.15 (the exponent was also computed via MLE [59]). In the inset we show the distribution of inverse intervals g(1/Θ), which is nearly uniform, since it is related to G(Θ) by g(1/Θ) = G(Θ)Θ 2 ∝ Θ 2−γ . The uniformity of g(1/Θ) will play an important role further on in the modeling of the pressure fluctuations PSD.

Connection between the burst size and time distributions
We consider now the link between the burst size distribution N (n) shown in Fig. 3 and the burst time distribution G(Θ) in Fig. 4. LetȦ = s/Θ denote the average growth rate of a burst of area s during the time interval Θ. This corresponds to an external perimeter growth [57,58], thereforeȦ ∝ ul e where u is a characteristic front speed (set by the Darcy law and the characteristic capillary pressure) and l e is the external perimeter, related to the linear size across a cluster l as l e ∝ l De . Since s ∝ l D , we have with β = 1 − D e /D. The distributions of s and Θ are linked by |G(Θ)dΘ| = |p(s)ds| and since the area s of a burst is proportional to its number of pores n (see Fig. 3), it follows that p(s) ∝ s −τ . Therefore, 6 Fluctuations in the measured pressure signal Next, we analyze the fluctuations in the pressure signal, following the pore invasion events. In Fig. 5, we show the typical pressure signature in a CIP experiment. The observed pressure pulses present a characteristic exponential relaxation. We also observe that a pulse can trigger others and even give rise to large avalanches with the invasion of several pores. A pulse can be divided into two phases: an initial fast drop in the capillary pressure p c and a slower exponential relaxation back to the pressure level ρgh set externally (see Fig. 1). The fast drop in p c occurs as the liquid is displaced (following the invasion of one or more pores) and subsequently redistributed to the surrounding menisci, causing a back-contraction of the interface [3,4,18]. The relaxation phase occurs as the liquid-air interface readjusts itself inside the where p w is the pressure in the wetting phase (liquid) just after the liquid-air interface, u is the average Darcy velocity of the flow in the porous network, µ = ρν is the liquid's dynamic viscosity, L i and k i with i = {1, 2, 3} are the length and permeability respectively of the porous network, filter and the tubing and S 1 and S 3 are the respective cross sections of the model and the tubing. Since the capillary pressure across the liquid-air interface is where is equivalent to an effective resistance to the flow. The volumetric flux in a pore is dV /dt = ua 2 /φ, where a is a characteristic pore length scale (for example the bead diameter) and φ is the porosity of the model. By introducing the concept of a capacitive volume κ = dV /dp c (used first in Ref. [3]), where dV is the liquid volume displaced from a pore throat in response to a change dp c in capillary pressure, we have Plugging this equation into Eq. (5), thus producing the exponential behavior seen in Fig. 5. C = p c (0)−ρgh < 0 is a constant associated to how much the capillary pressure decreases during the invasion of a set of pores before it starts to rise again. The characteristic time scale of the exponential decay is The invasion of one pore quite frequently triggers the invasion of others, in such a manner that before an exponential pulse decays completely, another one is seen in the pressure signal, see Fig. 5. This mechanism delays the complete relaxation of the pressure, effectively increasing the decay time from t c to t * ≥ t c . Indeed, if this relaxation-delaying mechanism was absent, the burst time distribution G(Θ) shown in Fig. 4 should be peaked around the value Θ = t c . Since we have shown that G(Θ) ∝ Θ −γ we expect the effective exponential decay time t * to follow the same distribution and, in particular, the effective decay rate λ = 1/t * should be uniformly distributed in an interval [λ min , λ max ] following the same distribution as 1/Θ (see inset of Fig. 4). λ max is related to the minimum decay time t * , i.e., λ max = 1/t c and we will consider λ min = 0 for convenience. Later on we will show that the distribution of decay rates has a crucial impact on the power spectrum of the pressure signal.

Pressure signal PSD
In Fig. 6 we show the power spectral density (PSD) associated to the pressure signal for the CIP experiments. The PSD S = S(f ) was computed for all experiments using the Welch method [64]. We have noticed the existence of a 1/f scaling regime (flicker/pink noise) for lower frequencies, followed by a crossover and a 1/f 2 scaling regime (brown noise) for intermediate frequencies.
For higher frequencies, another crossover follows and a region independent of f is seen (white noise associated with fluctuations in the pressure sensor and unimportant to our analysis). We see from Fig. 6 that the scaling properties of the power spectrum, in particular the occurrence of 1/f noise, seem to remain unchanged despite the changes in both sample dimensions and poresize distribution (the samples were rebuilt before each experiment, thus changing the pore-size distribution [56]). The 1/f regime is associated with events having frequency f < 10 −2 Hz, or alternatively, periods T > 100s. From Fig. 5, we see that this corresponds to the characteristic time intervals between the pressure pulses, thus indicating that they are associated with the presence of the 1/f scaling in the PSD.

Analytical modeling of the pressure signal and PSD scaling explanation
The non-trivial scaling of the CIP power spectral density can be explained by the following mathematical framework, which is an adaptation of an argument proposed in [65] to explain a similar 1/f to 1/f 2 transition in the very first reported observation of 1/f noise [47] (see also [66] and [67]). Apart from a nearly constant offset, the pressure signal can be modeled as a train of exponentially decaying pulses located at randomly distributed discrete times t j , where λ > 0 and A < 0 are initially taken to be constants (the characteristic decay rate and amplitude of the pulses) and H(t − t j ) is the Heaviside step function, i.e., H(t − t j ) = 0 if t < t j and H(t − t j ) = 1 if t ≥ t j . Let P λ (f ) be the Fourier transform of p λ (t). The PSD S λ (f ) is where r is the average rate of occurrence of pulses and the brackets are the expected value operator (since in practice one does not have access to an ensemble of measurements, we have employed Welch's method [64] to estimate the PSD, which is based on the concept of a periodogram [68]). The PSD shown in Eq. (11) is a Lorentzian curve which is approximately constant for lower frequencies (f ≪ λ/2π) and decays as 1/f 2 for higher frequencies (f ≫ λ/2π).
A model with a single constant decay rate λ cannot incorporate the 1/f region but, as previously argued, we expect λ to follow the uniform distribution ξ(λ) = 1/λ max in the interval [0, λ max ]. Taking this distribution into account and writing λ max = 2 π f t , we have (12) Eq. (12) has the asymptotic behavior thus presenting the 1/f to 1/f 2 transition observed in the experiments. The transition frequency f t in experiment CIP-1 is roughly f t = 1.5 · 10 −2 Hz (see Fig. 6). By using the constant A 2 r as a fitting parameter we can compare the measured PSD with the theoretical prediction in Eq. (12). Fig. 7 shows the resulting comparison produced using A 2 r = 1.5 P a 2 /s.
The dashed red vertical line marks the transition frequency f t . The analytical result reproduces the experimental findings very well, scaling as 1/f for f ≪ f t and as 1/f 2 for f ≫ f t . Indeed, this theory not only captures the 1/f and 1/f 2 domains but also fits the data well for the crossover region between these domains.
The transition frequency f t can be estimated using Eq. (9) and the resistance R from Eq. (6). As a first order approximation, let us consider only the contribution to R from the term R 1 relative to the resistance in the porous medium itself. Using µ = ρν = 5.1 · 10 −2 P a.s, L 1 = 0.27m, a = 10 −3 m, κ = 1.1·10 −12 m 3 /P a (from Ref. [3]), k 1 = 1.6·10 −9 m 2 and φ = 0.63 (both measured in a similar model in Ref. [5]), we find not far from the transition frequency f t = 1.5·10 −2 Hz shown in Fig. 7. The overestimation comes from the terms R 2 and R 3 in Eq. (6), ignored in the calculation above.
Finally, notice also the existence of a single isolated point in the very low frequency part of the PSD, falling far from the scaling region (extreme left for all experiments in Fig. 6). This point is not an outlier in the data: its existence signals the very slow positive drift of the pressure signal, which occurs since the capillary pressure has to increase to allow the invasion of narrower pores [3,56].

Comparison with a system driven under a CWR boundary condition
In order to test the effect of the boundary conditions in the PSD, we have run a controlled withdrawal rate (CWR) experiment using model (1). The resulting PSD is shown in the inset of Fig. 7. The PSD still presents an interesting scaling, but with different scaling regimes: 1/f 1.5 , for lower frequencies, and 1/f 3.5 , for intermediate frequencies.
The 1/f region is only observed for systems driven under the CIP boundary condition. The fact that the exponents for CWR differ from CIP is not surprising, since the pressure relaxation in that case is no longer exponential, but linear, see Ref. [3].
10 Connection between the measured pressure and the capillary pressure The pressure sensor measures the difference between the pressure in the air and the liquid at the outlet, i.e., p m = p nw − p out w . The measured signal is not exactly the capillary pressure p c = p nw −p w across the liquid-air interface, since p w = p out w given that viscous losses occur between the liquid-air interface and the outlet, thus generally making p w > p out w . Those losses occur in the porous medium itself and in the filter at the outlet of the model (numbers (1) and (2) in Fig. 1). The connection between p m and p c is p m = p c + u (R 1 + R 2 ), where R 1 and R 2 are the resistance terms from the porous network and the filter. Using Eqs. (7) and (8), we have Therefore, by comparing Eqs. (8) and (15), we see that p m differs from p c only in the amplitude of the pulses, but not in their characteristic exponential decay. Since our analysis depended only on the distribution of the decay rates, the differences between p m and p c are not crucial.

Further generalizations of the PSD analytical framework
One possible generalization of the model would be to consider a system with a distribution of amplitudes A instead of a single value (as we might expect from Fig. 5). In this case the scaling properties of the PSD would still be left unchanged but the constant A 2 in Eq. (12) and (13) would be replaced by the expected value of A 2 . Another possibility would be to consider a distribution for λ of the form ξ(λ) ∝ λ −δ . Here the 1/f 2 region is still left unchanged but the 1/f scaling is changed to 1/f (1+δ) [69]. As previously noted, the distribution of decaying rates λ is the crucial figure behind the 1/f scaling.

Conclusions
We have analyzed the burst dynamics from slow drainage experiments in porous media. We showed that this dynamics presents many features commonly associated to critical systems. Intermittent bursts of activity were observed over many time and length scales and a theoretical expression for their size distribution scaling, Eq. (1), was verified experimentally.
The pressure signal of the invasion presented an interesting PSD scaling, with a 1/f scaling region which further transitions to 1/f 2 in the case of the CIP boundary condition. We have employed an analytical framework [65] which satisfactorily reproduces the scaling properties of the PSD. The derivation of closed expressions relating the pressure signal PSD to properties of the porous medium and the fluids can lead to new techniques for indirectly probing such systems. For example, if one has access to the PSD only and not to the full pressure signal, the transition frequency f t can still be measured and information on the ratio k 1 /φ between the permeability and the porosity of the medium can be found via Eq. (14). If the PSD and f t are known, Eq. (12) can be fitted to measure the product A 2 r between the amplitudes and rate of occurrence of bursts.