Intermittent Dynamics of Slow Drainage Experiments in Porous Media: Characterization Under Different Boundary Conditions

The intermittent dynamics of slow drainage flows in a porous medium is studied experimentally. This kind of two-phase flow is characterized by a rich burst activity and our setup allows us to characterize those bursts directly via images of the flow and pressure measurements. Two different boundary conditions were analyzed: controlled withdrawal rate (CWR) and controlled imposed pressure (CIP). We have characterized geometrical and statistical properties of the bursts from images and pressure measurements. We have shown that in spite of leading to similar final invasion patterns, some dynamical features of the invasion differ considerably between the CWR and CIP boundary conditions. In particular, their pressure signatures are very distinct, which then translates into very distinct features on the power spectrum density of the pressure signals. A fully integrable analytical framework is presented which successfully describes the scaling features of the power spectrum for the CIP case.

Another interesting aspect of this problem is the fact that the dynamics of fluids in porous media very frequently presents intermittent behavior [3,19,23,24,53], with long intervals of stagnation followed by short intervals of strong activity in an irregular aperiodic fashion [56][57][58]. Intermittent behavior has been observed in various physical, biological, and economical systems [20,[59][60][61][62][63][64][65][66][67][68][69][70][71]. The occurrence of intermittent phenomena in vastly different fields points to the fact that its origin is typically not connected to specific details of the systems. Instead, intermittent behavior is generally associated to a competition between an adaptive external force driving the dynamics and an internal random resistance against that force. In the case of flow in porous media, the external force could come for example from an applied pressure difference across the medium, while the internal resistance is caused by the narrower or broader porethroats, whose capillary pressure thresholds need to be overcome to give continuation to the invasion [3,19,53].
In the present work we will show experimental results on the morphology and dynamics of invasion bursts during twophase flow in porous media. The flows studied are slow enough to be in the so-called capillary regime, in which capillary forces are typically much stronger than viscous ones [3,5]. We will employ artificially developed systems driven by two kinds of boundary conditions: controlled withdrawal rate (CWR) and controlled imposed pressure (CIP). Our experimental setup allows us to characterize the dynamics both by direct imaging of the flow and by local pressure measurements. For the CIP boundary condition, we present results related to the statistics of bursts, their morphology and orientation within the medium, as measured directly from the images of the flow. For the CWR boundary condition, we focus on the statistics of bursts defined over the pressure signal, following established literature conventions [72]. We will show that in spite of having similar overall morphology, the invasion dynamics is sensitive to the kind of boundary condition driving the flow. The pressure signature is clearly different and this difference is translated into the power spectral density (PSD) associated with the fluctuations in the measured pressure that follow the invasion events. It has been shown that for the systems driven by the CIP boundary condition, the pressure signal PSD presents an interesting 1/f scaling regime (pink noise) [23]. The observed PSD scaling is further described by employing an analytically integrable general mathematical framework, which, when reduced to a specific set of conditions, predicts both the 1/f scaling for lower frequencies, and a transition to 1/f 2 scaling for intermediate frequencies (brown noise), also observed experimentally. The general formulation of the analytical framework for the PSD has the potential to find applications in a much broader class of problems than the specific systems discussed here.

Description of the Experiment
The experimental setup employed is shown schematically in Figure 1A. We used a modified Hele-Shaw cell [73], filled with a monolayer of glass beads which forms the porous network. The beads have diameters a in the range 1.0 < a < 1.2 mm. The cell is made of one rigid rectangular Plexiglas plate (top plate) against which a layer of contact paper is attached with the sticky side facing away from the plate. The beads are randomly thrown onto the contact paper and get attached upon contact with the glue. After thoroughly filling the contact paper surface with beads, the system is shaken to ensure that the beads form a monolayer. A filter made of a spongeous material with small pores is placed close to the outlet of the model (the purpose of this filter will be soon explained). Silicon glue is used to define the in-plane boundaries of the model and another contact paper sheet is placed on top of it (with the sticky side facing the beads). The porous matrix is therefore bounded on top and bottom by the contact paper layers and on the sides by the silicon glue. The whole system is placed on the top of a pressure cushion, formed by a membrane filled with water connected to a reservoir placed at a fixed height of ca. 3 m. This cushion presses the beads (sandwiched in-between the contact paper sheets) against the upper Plexiglas plate and ensures that the beads won't move during the experiment, thus keeping the quasi-2D geometry of the medium despite the variations in bead diameter. In order to make the whole system more robust and avoid bending of the Plexiglas plate, an additional heavy glass plate is placed on top. Channels for the inlets and outlets had been previously milled in the Plexiglas plate and cuts were made in the upper contact paper sheet such that the liquid can be injected into and withdrawn from the porous network. Since the cell is placed horizontally, gravitational effects to the flow can be disregarded. From top to bottom the layers of the system are: glass plate-Plexiglas plate-contact paper sheet-porous matrix (glass beads)-contact paper sheet-pressure cushion-Plexiglas plate (see Figure 1A). Images of the flow are taken from the top by a digital camera (NIKON D7100) and pressure measurements are made at the outlet of the model with electronic pressure sensors (Honeywell 26PCAFG6G). The porous matrix was initially filled with a wetting viscous liquid composed of a mixture of glycerol (80% in weight) and water (20% in weight) whose kinematic viscosity, density and surface tension were measured to be, respectively, ν = 4.25 10 −5 m 2 /s, ρ = 1.205 g/cm 3 and γ = 0.064 N/m. Since the medium is initially completely wet by the liquid, the contact angle is always found to present low values (measured within the liquid phase), although its exact value varies during a pore invasion event due to dynamical effects.
In each experiment, air (non-wetting phase) invades the porous medium from the inlet channel, thus displacing the liquid (wetting phase) previously filling the pore spaces. Since in our experiments the non-wetting phase is the invading phase, we call them drainage experiments (in opposition to imbibition experiments in which the wetting phase invades the medium). Figure 1B shows a typical pore invasion event (imaged with a macro lens to enhance details). The invasion of a given pore happens when the capillary pressure p c , i.e., the pressure difference between the non-wetting and wetting phases, p c = p nw − p w , is high enough to overcome the capillary pressure threshold p t c associated with the geometry of the pore-throat giving access to that specific pore. Narrower pore-throats present higher values of capillary pressure threshold in accordance to the Young-Laplace equation [31], and are, therefore, harder to invade. During the dynamics, the capillary pressure builds up until the invasion of the widest pore-throat available to the airliquid interface happens (the one having the lowest capillary The liquid outlet on the right is connected to one of the modules shown in the dashed squares, depending on the choice of the boundary condition (controlled imposed pressure-CIP, or controlled withdrawal rate-CWR). In the CWR case, the outlet is just connected to the syringe in the pump, and in the CIP case, it is connected to the external reservoir placed above the moving platform. On the top we show a typical image of the experiment and a detailed zoomed in section showing the structure of the porous network, the invading air phase (white) and the defending liquid phase (blue). (B) A typical pore invasion event, imaged with a macro lens. Air invades the defending liquid through the widest pore-throat available to the liquid-air interface at that time (marked by the red arrow on the left image). After the invasion, a new set of pore-throats becomes available to the interface. Adapted from Moura et al. [23] under the Creative Commons CCBY license. pressure threshold). Once the pore connected to that pore-throat has been invaded, a new set of pore-throats becomes available for invasion and if one (or more) of these throats happens to have sufficiently low values of capillary pressure threshold p t c , they can also be invaded without further increasing the capillary pressure, giving rise to the characteristic burst dynamics we study [3,46,74]. These invasion events are followed by characteristic pulses in the measured pressure signal, having a typical relaxation signature which is nearly linear in time, for the CWR boundary condition, or exponential, for the CIP boundary condition, as we shall see in details in section 3.3.
The placement of the filter at the outlet boundary allows for the invasion to continue within the model even after breakthrough is achieved, i.e., even after the air phase first percolates through the model, forming a sample spanning cluster. The invasion continues until the moment in which all pores connected to the filter are invaded and the non-wetting fluid saturation reaches its final residual value. We recall that the saturation S of a given phase is understood as the ratio between the volume occupied by that phase in the porous medium and the total pore volume. The placement of the filter allows for a better statistics of the bursts in the vicinity of the outlet region, since the burst count in that area would be much reduced in the case without the filter in which the invasion stops once breakthrough is achieved. Next we move on to the description of the different boundary conditions driving the flow.

Description of the Boundary Conditions
Two distinct kinds of boundary conditions were used to drive the dynamics of the flow. We will call them controlled withdrawal rate (CWR) and controlled imposed pressure (CIP) boundary conditions. In this section we will describe them, starting from the latter.

Controlled Imposed Pressure (CIP)
The controlled imposed pressure (CIP) boundary condition was designed in order to drive the fluid invasion in a quasi-static manner such that the system would evolve from one equilibrium configuration to another in a slow, controlled fashion. In order to achieve that, the pressure difference across the model had to be dynamically controlled in such a way that it would be just enough to drive the invasion, but not exceedingly high. Had the imposed pressure difference been too high, the invasion would happen fast, and unwanted viscous and inertial effects would be present [75], removing the system from the capillary regime that we propose to study.
We have designed a mechanism composed of an external liquid reservoir (whose detailed description we postpone to Appendix 1 in Supplementary Material) and a logical feedback loop in order to control the pressure difference across the model. The outlet of the system is connected to the liquid reservoir by means of a continuous column of liquid. This reservoir is placed on the top of a moving platform whose height can be electronically set using a step-motor connected to a computer. By controlling the height of the reservoir, one can effectively control the pressure in the liquid phase (the pressure in the air phase is constant and equal to atmospheric pressure since the inlets are kept open in all experiments). The moving platform, liquid reservoir, and an extra syringe pump (whose detailed functionality is described in Appendix 1 in Supplementary Material) form what we call the CIP module, shown in the bottom dashed square on the right of Figure 1A.
Once the hydrostatic pressure caused by the height difference h (see Figure 1A) is enough to overcome the minimum value of capillary pressure threshold associated with the line of poresthroats at the inlet, i.e., once ρgh > min(p t c ) where g is the gravitational acceleration, the dynamics is triggered and the pore associated with the throat having the lowest capillary pressure threshold p t c is invaded. The liquid-air interface has now access to a set of new pore-throats, and if their values of p t c are low enough (i.e., if they are wide enough) the invasion progresses until the moment in which all pore-throats available to the interface have capillary pressure thresholds higher than the imposed capillary pressure, i.e., until the condition p t c > ρgh holds true for all available pore-throats. At this moment we say that the system has reached an equilibrium configuration and one has to increase the imposed capillary pressure by lowering again the position of the liquid reservoir (thus increasing the height difference h and the imposed capillary pressure p c = ρgh) up to the moment in which the invasion starts again. We define a burst in the CIP driven system as any connected set of pores and porethroats invaded between two time instants in which the imposed capillary pressure had to be increased.
The invasion dynamics is characterized by intervals of inactivity punctuated by localized burst events. This stick-slip process repeats itself throughout the experiment and it would be unpractical (and also very prone to errors) to keep a manual control of the imposed capillary pressure (the external reservoir height). In order to tackle this issue, we have devised an imagebased logical feedback loop to handle the control. Figure 2 shows a schematic of this feedback loop. The images recorded from the flow are thresholded to isolate the area A of the air phase. If this area is growing, as in the interval I-II in the figure, in which the area goes from A = A 0 to A = A 0 + dA, more pores are being invaded and we keep the imposed capillary pressure at the current level, say p c = p 0 . If, on the other hand, we perceive that the area growth didn't happen for a long enough interval (a bit longer than 1 min in our case), we say the front has reached a stationary configuration and in order to keep driving the dynamics, we increase the capillary pressure to a value p c = p 0 + dp. This is done by lowering the liquid reservoir by a given value dh related to dp by dp = ρgdh.
We record images every 15 s during the whole fluid invasion, and the image analysis for the feedback loop is done every 5th image. That means we compare the area of the air phase in, say, the 10th image with that of the 5th in order to decide between the states moving/not moving the liquid reservoir. The amount dh by which the reservoir is lowered is one of our control parameters. During the initial phase of the experiment, the liquid level in the reservoir is at the same height as the porous medium and the system needs to build up a considerable pressure until the invasion starts, therefore, we have chosen a larger dh = dh max for this initial process and a smaller dh = dh min to be used after the invasion has started (after the invading phase area reaches a certain small threshold). The values used were dh max = 100 µm and dh min = 10 µm corresponding to respective increments in the imposed capillary pressure of approximately dp max = 1.2 Pa and dp min = 0.12 Pa. 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 specific value of dh.
One particular point of interest for driving the system under the CIP boundary condition is the fact that in this case the invasion falls into the general class of problems in which an interface moves under the influence of an external (adaptive) force through a medium with quenched disorder. These problems have received a great dose of interest from the scientific community and have been extensively studied theoretically and numerically. Special attention has been paid to the transitions that occur as the external force is increased up to a critical threshold and several scaling relations between the exponents characterizing these critical transitions have been proposed [76,77]. Our experimental setup thus allows us to directly test the applicability of some of those relations in porous media flows [23].

Controlled Withdrawal Rate
The controlled withdrawal rate (CWR) condition is relatively simpler to implement. Its experimental realization is done by replacing the CIP module in Figure 1A by a syringe pump attached to the model's outlet. There is no need for a feedback loop to control the experiments. In the experiments reported here, we used a constant slow withdrawal rate of q = 0.005 ml/min, thus assuring that the system remains in the capillary regime [5]. The associated capillary number is given by Ca = ρνq/ γ = 6.1 · 10 −7 , where = 1.1 · 10 −4 m 2 is the crosssection of the model. The CWR boundary condition is applied in the majority of two-phase flows studies present in the literature [3,19,49,52]. For the CIP we cannot compute the Capillary number using the same formula, since the flow rate q is not constant. Nevertheless, since the duration of both experiments is similar, we can assume their capillary numbers to be comparable.
When applying a syringe pump to drive the system, one important detail to be considered is the choice between FIGURE 2 | Feedback loop used to drive the CIP experiments. This figure shows a detail from the black and white thresholded image, where the air phase is white and the liquid phase together with the porous medium are black. From I to II, the area of the air phase has grown from A 0 to A 0 + dA, so the capillary pressure is kept constant at p 0 (no change in the liquid reservoir position). From II to III, the area of the air phase does not change, which means the capillary pressure must be increased from p 0 to p 0 + dp (by lowering the liquid reservoir and thus decreasing the pressure in the liquid phase). This analysis is done "on the fly" as the experiment is performed. Adapted from Moura et al. [74] under the Creative Commons CCBY license.
withdrawal of the wetting phase or injection of the non-wetting phase (in the case of drainage studied here). If the compressibility of both phases is negligible, as, for example, in the case of slow displacement of one liquid by another, these two options lead to similar results. Nevertheless, in the case in which one of the phases is compressible, it is preferable to have the syringe pump connected to the other phase, in order to avoid additional effects to the dynamics arising from compressibility (unless those effects are precisely the object of the study, see Jankov et al. [78]). Since air is compressible, we made the choice of having the syringe pump withdrawing the liquid instead.
Another important aspect of this boundary condition is the fact that the liquid displaced from one pore when it is invaded by air cannot immediately leave the system. This happens because the fixed withdrawal rate and typical flux rates during the bursts have very different time scales. In the time scale of a single pore invasion, the volume available for the fluid in the syringe is essentially unaltered. Since the displaced fluid has to go somewhere, it is redistributed to the liquid-air interface in the pore-throats in the vicinity of the invaded pore. This important interface back-contraction effect is responsible for the Haines jumps [53] in the capillary pressure following each pore invasion and has been successfully incorporated into some invasion percolation models to simulate two-phase flow in porous media [3,19].

Dynamics of the Air Saturation
To give an idea of the evolution of an experiment, we show in Figure 3A a sequence of snapshots of the flow (here driven by controlling the imposed pressure-CIP boundary condition). The current time (in hours) and imposed capillary pressure (in Pascals) are shown under each image. The time difference between the first and last images is about 82 h. The air inlet is on the left and the filter is the black stripe on the right. In Figure 3B, we show a spatiotemporal map of the full invasion. The color map gives the invasion time of each pore in the system. Notice that if the outlet filter were not present, the invasion would have stopped at the moment of breakthrough, when the air phase first percolates through the model. This corresponds to t ≈ 60 h in Figure 3B. Without the filter, the whole red region in the figure would not have been invaded.
The spatiotemporal map shown in Figure 3B is a very useful matrix: once this map is computed, one does not need to rely on a large collection of pictures/video any longer. In practice this brings the benefit of dramatically reducing the computational time for any further analysis.
Since we are controlling the imposed pressure (by controlling the level of the external reservoir, see Figure 1A), we can also define a pressure map for the flow, here shown in Figure 3C. The color code (shifted with respect to Figure 3B to aid visualization) indicates the imposed capillary pressure during the invasion of a given pore. Notice that for the system studied here, the majority of the invasion happens in a relatively narrow band of imposed pressures between 420 Pa and 480 Pa.
The box counting fractal dimension [79] of the invading structure was measured to be D = 1.76 ± 0.05 for CIP and D = 1.75 ± 0.05 for CWR, therefore, with respect to the morphology, one can say that both boundary conditions yield similar structures (for similar capillary numbers). Nevertheless, despite these similarities in morphology, there are some very clear differences in the dynamics of the different systems. In particular, the evolution in time of the air phase saturation happens very differently. For the CWR, since the withdrawal rate of the liquid is constant, externally set by the syringe pump, we expect the area growth of the air phase to happen linearly, in a steady controlled manner, following the same rate of the imposed withdrawing. This behavior is observed in Figure 4, where the phase saturation is determined directly from the images of the flow. The red curve (bottom) shows the data for the CWR case. We see that once the invasion starts, it progresses at the constant rate imposed by the syringe pump. Conversely, for the CIP boundary condition, the air saturation evolution happens in a FIGURE 3 | (A) Images from the drainage process driven by the CIP boundary condition. The time t of each image is shown in hours and the instantaneous capillary pressure p c in Pascal. The porous medium is initially saturated with the liquid phase (blue) which is then gradually displaced by the invading air phase (white). On (B) we show a spatiotemporal map of the invasion, the color code indicating the invasion time of a given pore. On (C) we show the imposed pressure map, in which the color code denotes the imposed pressure (given by the height of the external reservoir), at which the invasion of a given pore occurred. much more intermittent fashion, in which intervals of stagnation are followed by periods of sudden growth, as shown in the blue curve (top) in Figure 4 and more clearly in the zoomed-in inset. These differences in the dynamics of invasion will be reflected later in the power spectral density of the pressure signal, to be analyzed in section 3.5.

Geometrical Characterization of the CIP Bursts
Our setup allows us to directly visualize the invasion bursts from the images of the flow. Let us consider the experiment shown in In an ideal scenario, with infinite precision, once the system reaches an equilibrium state, i.e., one at which the imposed pressure difference is not enough to invade any of the pore-throats available to the liquid-air interface, we would like to be able to increase the pressure difference ever so slightly to trigger the invasion through one and only one of the pore-throats (the widest among the ones available to the liquid-air interface). A burst could then be naturally defined as the region that is invaded following the invasion of this particular pore-throat until the interface reaches a configuration in which again all FIGURE 4 | Temporal growth of the air saturation measured from the area of the invading air phase (same model shown in Figure 3). Time and saturation are normalized by their final values, t max and S max , respectively. The blue curve (top) shows the data for the controlled imposed pressure (CIP) boundary condition, whereas the red one (bottom) is for the controlled withdrawal rate (CWR). As expected, the later shows a rather constant growth rate, set from the external syringe pump. For the CIP case, we see a much different scenario, with intervals of relative stagnation being followed by periods of sudden area growth in an intermittent manner. The zoomed-in insets emphasize this difference.
pore-throats available are too narrow to be invaded and the imposed capillary pressure has to be increased once again. This burst's size would be given by the number of invaded pores and the burst's duration could be identified as the time interval between the two instants at which the imposed pressure was increased (immediately before and after the burst). With our current setup, with pressure increments of the order of dp min = 0.12 Pa, we still don't get to that idealized level of accuracy and very frequently, after increasing the pressure, the invasion is triggered from a small number of pore-throats (between 3 and 4 on average). We therefore consider the connected invasion starting from each of these pore-throats as separate bursts, all of which are associated to the same time interval . Our definition of bursts for the CIP experiments is then the connected set of pores that have been invaded from a given pore-throat until the imposed pressure needs to be again increased. Notice that this definition makes explicit use of the fact that the bursts happen at constant imposed pressure, therefore, it is only suitable for the CIP case. Figure 5 shows the resulting bursts. In the top part of the figure, they are colored by their area, normalized by a characteristic pore area (100 pixels ≈ 0.3 mm 2 in these images). In the bottom part of the figure, the same bursts are shown with random coloring, in order to emphasize their individual shapes. Notice that the vast blue areas in the top image are indeed formed by a large number of small bursts (see the detail in the images). Only bursts having their centroids lying in the central 90% of the model's length are considered, in order to avoid possible boundary effects which have been observed to happen in the vicinity of the inlet and outlet of the flow [74].
The experimental assessment of individual bursts provided in Figure 5 yields a valuable amount of information about the dynamics. It has been shown [23] that the distribution N(n) of bursts of size n (where n is the estimated number of pores, computed by normalizing the burst area by a typical pore area) follows the scaling N(n) ∝ n −τ , where τ = 1.37 ± 0.05. The burst dynamics is therefore self-similar, having no characteristic intrinsic length scale [80]. The exponent τ measured is consistent with the theoretical value τ = 1.30 ± 0.05 predicted by numerical simulations (invasion percolation) and analytical results [23,45,76,81].
From Figure 5 one cannot directly notice any tendency of concentration of larger bursts toward one of the ends of the model. Indeed, the amount of bursts and their average sizes seem to be approximately homogeneously distributed along the system's length. In order to validate this statement quantitatively, we have analyzed the average burst sizen s (burst area divided by typical pore area ≈ 3 mm.) and the total number of bursts ξ s within sliding subwindows having a fixed length of 2 cm and width spanning the whole width of the model, see Figure 5. For each subwindow, we have computed the number of bursts ξ s whose centroid lie inside it and their corresponding average bursts sizen s . The results are shown in Figure 6A. The correlation observed between consecutive measurements arises from the fact that neighboring subwindows overlap: the centers of consecutive subwindows are separated by 0.2 cm, whereas the subwindow length is 2 cm, thus yielding 90% overlap between neighboring subwindows. From the analysis of this figure, we cannot observe any systematic trend for the average size or number of bursts: when seen separately, neither of the average burst size nor the number of bursts seem to be systematically increasing or decreasing with the subwindow position. Nevertheless, when comparing one against the other, we can observe a clear anti-correlation of the data. This anticorrelation is confirmed by measuring the Pearson's correlation coefficient, given by wheren s i and ξ s i correspond respectively to the values of the average burst size and number of bursts inside the ith subwindow (i = 1, 2 . . . N), . denotes averaging and σn s and σ ξ s are the standard deviations of the respective quantities. Since the Pearson's coefficient is negative, the average burst sizen s and number of bursts ξ s measured in a sliding subwindow are anti-correlated. The origin of this anti-correlation is easily understood: if the system is in a steady regime, the total invaded area in a subwindow is approximately constant (apart from statistical fluctuations arising from the randomness of the porous medium and the underlying dynamics [74]). The products of the two quantities, average burst sizen s and the number of bursts ξ s inside a subwindow, gives the total area of bursts having centroid inside that subwindow. This is not exactly equal to the total invaded area A s inside the subwindow, since the boundaries of the bursts (particularly the larger ones) could cross the imaginary subwindow boundary, but it gives an approximation to it. Therefore, one can expect the product between the average burst size and the number of bursts inside As mentioned previously, the dynamics of our drainage experiments happen in the capillary regime [3,5], i.e., it is governed primarily by capillary forces. Since the pore-size distribution is nearly isotropic, due to the construction routine employed to make the models in which the glass beads are randomly placed, the capillary forces involved should also be isotropically distributed. We expect this property to be translated into the dynamics, generating bursts with no privileged direction within the model. One can have a qualitative picture of this statement directly from Figure 5 and a more quantitative result can be obtained by producing a histogram of the bursts orientation angle θ . This angle is defined as the one between the x-axis (inlet-outlet direction, see Figure 5) and the major axis of an ellipse having the same second-moments as that particular burst. Figure 6B shows the resulting histogram. We can see that, apart from fluctuations caused by the small size of the dataset, there isn't any particularly privileged orientation for the bursts. Indeed, such a distribution also indicates that viscous effects are much reduced in the experiment: if the dynamics were driven by viscous forces, one would observe bursts much more elongated in the inlet-outlet direction, following the average gradient of pressure in that direction. A distribution of burst angles would not be isotropic, but show a peak toward smaller angles and a corresponding reduction for larger angles. Up to this moment, we have concentrated our analysis mostly on quantities that could be obtained from the images of the flow. In the coming sections, we will shift our focus to the analysis of the pressure measurements, for both boundary conditions. We begin by briefly analyzing the typical signature of the measured pressure signal during the invasion.

Typical Signature of the Pressure Measurements
As anticipated, in spite of having nearly indistinguishable largescale morphology, the dynamical features of the invasion differ considerably between the CIP and CWR boundary conditions. This difference is clearly observable in the evolution of the measured pressure signal, as shown in Figure 7. Let us focus initially on the CWR signal in Figure 7A. Clear pulses in the capillary pressure corresponding to the pore-invasion events are observed (also termed Haines jumps [53]). Each pulse is composed of one drop phase and one relaxation phase (see arrows in the figure). Morrow [82] has used the terms rheon and ison to describe drop and relaxation phases respectively. The  Figure 5). The sliding subwindow spans the whole width of the system and has a fixed length of 2 cm. The spatial distribution of bursts follows no systematic trend along the model. The average burst sizen s is computed by dividing the burst area by typical pore area ≈ 3 mm. (B) Histogram showing the random orientation of the bursts shown in Figure 5. The thicker green line shows the mean of the distribution and the thinner red lines show its standard deviation.
pressure drops correspond to the moments in which one or more pores are invaded and the liquid volume that previously filled them is redistributed to the menisci in the surrounding porethroats (see Figure 1B). The slower relaxation part corresponds to the moments in which the menisci of the air-liquid interface slowly advance inside the pore-throat, as a response to the constant withdrawal rate imposed by the external syringe pump. The capillary pressure slowly increases, as the menisci reach narrower parts of the pore-throat up to the moment when the pressure threshold associated to the widest pore-throat is reached, the interface then becomes unstable and a new poreinvasion occurs, giving rise to another pressure drop [3].
The connection between the capillary pressure and the volume of liquid redistributed to (drop phase) or withdrawn from the interface (relaxation phase) was explored in Måløy et al. [3] by introducing the concept of a capacitive volume κ such that dV p = κdp c , where dV p is the liquid volume displaced from an interface pore throat in response to a change dp c in capillary pressure. The volume dV p can be connected to the total volume dV redistributed to or withdrawn from the whole interface as dV p = dV/n f , where n f is the number of pore-throats belonging to the liquid-air interface. By putting these pieces together we get Considering that the process occurs in a time dt, we have that the variation dp c /dt is connected to the total flux q i = dV/dt by where q i = q R for the relaxation phase and q i = q D for the drop phase. During the relaxation phase the flux q R withdrawn from the interface, matches the externally imposed withdrawal flux (q R = q = 0.005ml/min in our case). For the drop phase the flux q D redistributed to the interface can be estimated as q D = n f a 2 v b , where v b is a characteristic back-contraction velocity of the interface. Like the Darcy velocity [1], v b is expected to vary inversely with the viscosity of the fluid. So for experiments conducted with a less viscous fluid, one would expect a higher v b which leads to a higher redistribution flux q D during the drop phase and consequentially a steeper pressure drop in accordance with Equation (4). Indeed this can be seen when comparing the drop phase of our experiment in Figure 7A (where a highviscosity liquid was used) against the pressure signal of an experiment in which water is used, like in Figure 1D in Måløy et al. [3] for example. Notice that in the figure in Måløy et al. [3], the vertical axis denotes the pressure in the liquid phase, while in our case we show the capillary pressure instead, i.e., the difference between the pressure in the air and liquid phase. So the downwards drops in Figure 7A are to be compared against the steeper upward sections of the signal in Måløy et al. [3]. While q R is an externally imposed parameter, q D is to be understood as an intrinsic parameter, dependent on the viscosity of the wetting phase.
Notice that the inverse dependence of the capillary pressure drop with the number of pore-throats in the system shown in Equation (3) implies that the existence of the pressure drops (Haines jumps) is a direct consequence of the finite size of the system. For a large enough system (n f − → ∞), the pressure jumps cease to exist [3].
If one zooms further in one of the relaxation phases, one can observe (left inset in Figure 7A) some low amplitude fluctuations in the signal, corresponding to small adjustments of the menisci inside the pore-throats that happen before the invasion of another pore. Such adjustments are not easily observable from our images, since they happen at very short time and length scales. Nevertheless, they are bound to occur in any case of slow flow in a porous medium (both in drainage and imbibition). Zooming even further in on the signal allows one to see the high frequency random fluctuations (right inset in Figure 7A) which are expected from any electronic measuring device, in our case the pressure sensor employed. These fluctuations set the error in the pressure measurements to be on the order of δp = 0.5 Pa.
For the CIP boundary condition, shown in Figure 7B, we see that the pressure signal presents an entirely different signature. Most of the time the system is in a stagnation state, in which the capillary pressure, set by the external reservoir position, is (very slowly) increased. This corresponds to the nearly horizontal level approximately at 405 Pa in Figure 7B. Those large intervals of stagnation are punctuated by sudden activity events, seen as the nearly vertical drops followed by the exponential relaxation phases in the figure. Energy is slowly fed into the system during stagnation periods between pore invasion events and is rapidly released (and further dissipated by viscosity) during the fast invasion events. Differently from the CWR, in the CIP case the pressure pulses present a nearly exponential relaxation behavior. One pulse can eventually trigger others or even give rise to a large avalanche of multiple pressure pulses, as seen in the inset of Figure 7B.
In the CWR system, the total available liquid volume in the experiment is composed by the volumes of (a) the liquid phase in the porous medium, (b) the volume inside the tubing connected to the outlet and (c) the (slowly changing) volume of the syringe pump. Once a pore is invaded by air, the liquid volume previously filling that pore has to be redistributed to the liquid-air interface in the pore-throats in the vicinity of that pore because the time scale needed for the pump to buildup that extra volume in the syringe outside the medium is much larger than the typical time scale of a burst. On the other hand, when the system is driven under the CIP boundary condition, after a pore invasion, the interface back-contraction is much reduced, since the liquid previously filling the pore can freely flow outside of the system and into the pressure reservoir (see Figure 1A). This difference in boundary condition is ultimately responsible for the nearly exponential relaxation behavior shown in the CIP case. A full derivation of this behavior was done in Moura et al. [23] and we include it here in Appendix 2 in Supplementary Material.
The differences in the pressure signatures is further reflected on the signal's power spectral density, as we shall see in section 3.5. In the following section we will perform an analysis of time-directed avalanches in the CWR pressure signal.

Avalanches in the CWR Pressure Signal
In this section we will give a statistical characterization of the pressure signal for the CWR experiments by considering the concept of time-directed avalanches. This concept was previously used in numerical simulations by Roux and Guyon [72] and Maslov [83] and in both simulations and experiments by Aker et al. [84] and Biswas et al. [85]. Figure 8A shows the full pressure signal and a zoomed in section is seen in Figure 8B. The forward time-directed avalanches are presented as the green horizontal lines, black dots denote the set of local maxima of the signal (starting points for the avalanches) and the red crosses correspond to the local minima of the signal. Each forward avalanche corresponds to the time interval T f between the moment at which a pressure drop starts and the first moment at which the pressure reaches a level equal to that from the starting point (in the future with respect to the starting time). Conversely, a backward avalanche can be defined as the time interval T b between a local maximum (again one of the black dots) and the last moment at which the pressure reached a level equal to that from the starting point (in the past with respect to the starting time). Backward avalanches can be computed simply as forward avalanches in the time-reversed pressure signal. We have included in Figure 8B one arrowed green line to denote a forward avalanche of size T f and one arrowed purple line to denote a backward avalanche of size T b .
Notice that the definition of time-directed avalanches employed here has a hierarchical structure: a particular avalanche can encompass a number n of others. Each avalanche defines a "pressure valley, " the region of the pressure signal between the beginning and end of the avalanche. Each valley is composed of either one (in the case of a singular avalanche) or a series of pressure drops (compound avalanche) each drop corresponding to a pressure difference p D i (difference between the pressure level of a black point and the following red point in Figure 8B). The forward valley size χ f is defined as the sum of all the pressure drops within the valley, i.e., Backward valley sizes can be defined in an analogous manner.
Frontiers in Physics | www.frontiersin.org For each avalanche of size T f , we can then measure the respective valley size χ f . In Figure 9 we show a scatter plot between the two quantities for all forward avalanches (singular and compound) in the experiment. We observe that there is a clear correlation between the quantities, i.e., χ f ∝ T f . The Pearson's correlation coefficient was measured to be ρ = 0.995. This correlation is not accidental and can be easily understood. Let us assume that the derivatives of the pressure signal in the drop and relaxation phases are approximately constant, given by Equation (4). Notice that during a given avalanche, the pressure goes down a number n of times and the total time spent in these drop phases is given simply by T D f = p D i / dp c /dt . A similar expression can be derived for the total time spent on the relaxation phases for a given avalanche, i.e., T R f = p R i / dp c /dt , where p R i denote the increments in pressure. Notice that the time derivative of the pressure signal is different for each phase, in accordance to Equation (4). Notice also that although the consecutive pressure differences associated with the drop and relaxation phases are not necessarily equal, the sum of all the pressure drops during an avalanche equals the sum of all pressure increments, i.e., p R i = p D i = χ f . The reason for that is the fact that subtracting all the increments from all the drops in pressure equals zero for any avalanche, since the end point of an avalanche has, by definition, the same pressure as the starting point (same height in Figure 8B). To put it differently, the cumulative lowering in capillary pressure equals the cumulative increments for any given avalanche, since the pressure at the end equals that of the beginning of the avalanche.
Using again Equation (4), we have that the total time spent in the drop phases for a given avalanche equals where q D and q R are respectively the volumetric fluid flux given to the interface during the drop phase and withdrawn from it during the relaxation phase. Notice that when using Equation (4) to obtain Equation (6), we have assumed the linearity relation between the pressure drop and the fluid volume given to (drop phase) or withdrawn from the interface (relaxation phase), as expressed in Equation (3). Equation (6) shows a neat connection between statistical quantities defined over the pressure signal (the avalanche size T f and valley size χ f ), the current state of the invasion (via the size of the invasion front n f ), an externally imposed control parameter (the flow rate q R which here matches the external syringe flow rate q) and intrinsic properties of the system (the capacitive volume κ and the flow rate q D , the later depending inversely on the viscosity of the wetting phase).
Next we analyze the size distribution of forward avalanches N all (T f ). For convenience, we will use the same notation N all (X) to denote the distribution of any quantity X. The subscript indicates that one counts all avalanches (singular and compound). The distribution N all (T f ) is shown in Figure 10A. We have measured the scaling N all (T f ) ∝ T −τ all f , with τ all = 1.67 ± 0.05 over times ranging from 100s (the typical duration of a short singular avalanche in our experiment, see Figure 8B), to around 20000s corresponding to a large compound avalanche. In the inset of the figure we present a scatter plot showing the correlation between the avalanche size and the number of sub-avalanches within it. As anticipated, larger avalanches are typically composed of a larger number of sub-avalanches, as we see in the plot. The measured Pearson's correlation coefficient was ρ = 0.998.
The scaling exponent for the distribution N all (T f ) measured here differs from value τ all = 1.9±0.1 reported in Aker et al. [84]. It also differs from the superuniversal theoretical prediction τ all = 2 for invasion percolation and a number of other growth models presented in Maslov [83]. There are a number of factors behind this difference. First, the standard invasion percolation algorithm does not include effects from viscous dissipation and the back-contraction of the liquid-air interface after each pore invasion, present in the experimental case. Additionally, the pore structure of our models is not entirely uncorrelated. Some degree of correlation must exist due to the technique employed in the model construction. As previously explained, the setup is constructed by letting the glass beads fall onto the sticky side of a contact paper sheet. If the glass beads were mathematical points without dimension, randomly placing them in the contact paper should ideally generate an uncorrelated structure, but this is not the case when the beads have a finite size (1mm diameter in our case) and cannot physically overlap. In a recent study, Biswas et al. [85] have shown that indeed very little correlation in the pore structure can lead to changes in the critical exponents associated with the invasion process. The authors have shown in particular that the forward avalanches exponent does not need to assume the superuniversal value τ all = 2. Indeed, using numerical simulations Biswas et al. have shown that the forward avalanche exponent assumes the value τ all = 1.66 ± 0.04 for correlated porous samples, which is in direct agreement with our measurements presented in Figure 10A. When the simulations were performed in an uncorrelated porous medium, the authors have measured the exponent τ all = 1.99 ± 0.05, in agreement with the theoretical predictions in Maslov [83].
As an additional check of our results, we have measured the valley sizes distribution N all (χ f ). Since the quantities T f and χ f are clearly correlated (as seen in Figure 9), the distribution N all (χ f ) should present the same scaling behavior as N all (T f ).
Indeed this is what we observe, see Figure 10B where the black line denotes the scaling χ −τ all f , with τ all = 1.67 ± 0.05. In the inset we show the scatter plot between the valley size χ f and the number n of sub-avalanches inside it. Not surprisingly, we find that χ f ∝ n since it had previously been shown that T f ∝ n and χ f ∝ T f (Figures 9, 10A).
There's an important technical point about the computation of the avalanches employed here. As previously seen, it requires finding the local maxima and minima in the pressure signal as an initial step. In a smooth signal, this can be done easily by simply finding the points where the first derivative of the time series changes from positive to negative for the local maxima and vice-versa for the minima. But in a natural noisy signal, this approach is problematic due to high frequency fluctuations. One possible way to handle the issue is by employing a low-pass filter to remove the high frequency noise. This was our initial approach but it led to biases in the statistics, particularly in the region of small valley sizes. We decided to avoid the filter completely and perform the analysis in the raw signal employing a different technique to find the local extrema. We have employed Matlab's findpeaks routine to find all maxima in the noisy data and sort them out with respect to the peak prominence (a measure of the distance between the top of a particular peak and the surrounding signal). This can be used to remove most of the peaks mistakenly assigned to the fluctuations because they are very small. Next, we have excluded from the data all peaks that were followed by a pressure drop smaller than a certain small threshold, which we have arbitrarily chosen to be p D = 4 Pa. This threshold is chosen arbitrarily, but we have checked that all exponents of the distributions reported here did not vary significantly for any value in the interval 2 Pa < p D < 6 Pa.
For completeness we also include the distribution of backward avalanches T b and valley sizes χ b . These are found in Figures 11A,B, respectively. The insets show the correlations between T b and n in (a) and χ b and n in (b), where n denotes again the number of sub-avalanches. Interestingly the distributions of T b and χ b seem to present a crossover behavior, where the small avalanches scaling show an exponent that is smaller than the one for the large avalanches. We have measured the scaling exponent τ all b = 1.55 ± 0.05 for avalanches larger than a given crossover scale and τ all b = 0.90 ± 0.08 for avalanches smaller than this scale. The existence of a very similar crossover behavior had been reported in simulations performed in Biswas et al. [85], where the authors explained that small avalanches are influenced by effects coming from the local correlations in the pore structure. Larger avalanches on the other hand, being averaged over many pores, do not feel such local correlations. Our experiments are in direct agreement with the exponents reported in Biswas et al. [85], where it was measured τ all b = 1.59 ± 0.06 for the larger avalanches and τ all b = 0.93 ± 0.03 for the smaller ones.
The distribution of starting pressures (black dots in Figure 8) for the avalanches can also be directly obtained from our data. It is shown in Figure 12. Notice that, the majority of the avalanches have starting points in the interval between 420 Pa and 480 Pa. This is in line with what we had previously observed in the CIP case ( Figure 3C). The measurement in Figure 12 is also similar to what was found for an experiment performed in a similar model with different fluids in Furuberg et al. [19].  Figure 8).
Next we move on to the analysis of the power spectral density of the pressure signals. The analysis presented here includes and extends the work initiated in Moura et al. [24].

Power Spectral Density of the Pressure Signal
The different temporal signatures for the CWR and CIP pressure bursts are translated into the power spectral density (PSD) of the signal. Figure 13 shows the PSD S = S(f ) for the CWR (red dots) and CIP (green crosses) boundary conditions, the later being shifted upwards 4 orders of magnitude to aid visualization. Both PSDs present a scaling behavior S ∝ f −α , with two apparent crossover regions. We show in Figure 13 a set of guide-to-eye lines with the respective slopes for each scaling regime. Let us focus initially on the CWR case. In order to analyze each of these scaling regions separately, let us look again at the typical signature of the pressure signal for the CWR experiment, shown in Figure 7A. On the larger scale, we see the characteristic pulses in pressure, each of which is caused by the geometrical invasion of one or more pores by the air phase. Those pulses happen at intervals that match the period range of the first scaling regime in the PSD (α = 1.5 for frequency f in the approximate interval [0.001, 0.005 Hz], with corresponding periods T in the range [200, 1,000 s]), which indicate that their presence is responsible for this scaling regime. If we zoom into that scale, we see (inset on the left in Figure 7A) that there are some oscillations of the pressure signal happening in the phase of pressure build up that precedes the bursts. These happen with typical periods in the range corresponding to the second scaling regime (α = 3.5 for f in [0.007, 0.018 Hz] and corresponding T in [50, 140 s]). We believe those pressure oscillations are related to tiny adjustments of the air-liquid front, happening in between burst events while the front is in a stable position, as previously stated. Such adjustments are not easily observable from our images, since they happen at very short time and length scales. Finally, if we zoom in even further on the signal, we see (inset on the right in Figure 7A) the typical high frequency random fluctuations expected from any electronic measuring device. In the case of our pressure sensor, they are on the order of 0.5 Pa and, as expected, FIGURE 13 | Power spectral density of the pressure signal for CWR (red dots) and CIP (green crosses) boundary conditions, the later being shifted upwards 4 orders of magnitude. Both systems present scaling behavior for the PSD of the form S ∝ f −α , with different exponents α shown by the guide-to-eye lines. In both cases, the spectra for high frequencies present the typical flat profile characteristic of white noise from the experimental measuring device. Notice that for the CIP case, a region consistent with 1/f scaling (pink noise) for the low-frequency events is observed, followed by a 1/f 2 (brown noise) region for intermediate frequencies.
lead to the high frequency white noise region of the PSD, the plateau on the right of Figure 13 (α = 0 for f in [0.5, 1.76 Hz] and corresponding T in [0.5, 2 s]). The highest frequency in the graph is the Nyquist frequency, f N = 1.76 Hz, which corresponds to half of the acquisition frequency employed in that particular CWR experiment. Notice also that there is a single isolated point in the very low frequency part of the PSD (extreme left in Figure 13) that falls far from the scaling region. This point is present in all experiments and is not an outlier in the data: its existence signals the very slow positive drift in the pressure signal, which occurs since the capillary pressure has to increase (albeit very slowly after the front reaches a propagating regime [74]) in order to give access to the invasion of narrower pores as the front progresses inside the medium.
Next, let us move to the analysis of the pressure signal associated to the CIP boundary condition. As previously mentioned, the CIP pressure pulses look very different from the ones associated to the CWR case. They are composed of a sudden drop followed by a nearly exponential relaxation behavior, which is rather different from the signature shown for the CWR case (see Figure 7). This difference is reflected in the power spectral density for the CIP case, shown in green (crosses) in Figure 13. We notice that once again the PSD presents the scaling S ∝ f −α with three characteristic scaling regimes, but this time with exponents α = 1 (pink noise) for the lower frequencies, α = 2 (brown noise) for intermediate frequencies and again α = 0 (white noise) for the higher frequencies. As earlier, the lines shown in the plot are to be understood only as guide-toeye lines and not as accurate measurements of the exponents (which would require data extending for a longer range of time/frequency scales). The scaling region related to the bursts now has exponent α = 1 (1/f noise, also known as pink noise). The appearance of 1/f noise and the following crossing over to 1/f 2 scaling is an interesting observation, which has been noticed in a number of other systems, including in the very first reported observation of 1/f noise [86] (flicker noise in vacuum tubes).
In the following section we present a fully integrable mathematical framework that can be used to explain the unusual scaling behavior observed for the CIP boundary condition. This framework is an generalization of the one presented in Moura et al. [23]. A similar framework for the CWR case is still lacking.
3.6. Analytical Model for the 1/f to 1/f 2 Transition in the CIP Power Spectral Density As noted in Figures 7, 13, the pressure signal in the CIP case presents bursts with an exponential relaxation behavior and yields an interesting non-trivial PSD presenting a 1/f to 1/f 2 transition (pink noise to brown noise transition). In this section we propose an analytical framework to explain this unusual PSD scaling. This framework is an adaptation of an argument proposed in Ziel [87]. A similar argument was also presented by Bernamont [88] and a review can be found in Milotti [89].
Motivated by Figure 7, let us consider that, apart from a nearly constant offset, the pressure signal can be modeled as a train of exponentially decaying pulses located at a set of randomly distributed discrete times t j where λ > 0 and A < 0 are initially assumed to be constants (the characteristic decay rate and amplitude of the pulses) and H(t − t j ) is the Heaviside step function [90], i.e., H(t − t j ) = 0 if t < t j and H(t − t j ) = 1 if t ≥ t j . The characteristic decaying time of the pulses is given by t c = 1/λ. Let P λ (f ) be the Fourier transform of p λ (t), The λ-dependent Power Spectral Density S λ (f ) can be computed as where r is the average rate of occurrence of pulses (average number of pulses per unit time) and the brackets correspond to the expected value operator. In the practical calculations that will follow, since we do not have access to an ensemble of realizations of the signal, we will use the Welch's method [91] to estimate the PSD, which is based on the concept of a periodogram [92]. The PSD shown in Equation (13) of a train of exponential pulses of equal amplitudes A and decaying rates λ is therefore a Lorentzian curve, given by Equation (13). This function is approximately constant for lower frequencies (2πf ≪ λ) and decays as 1/f 2 for higher frequencies (2πf ≫ λ). Therefore, the model of a train of equally decaying exponential pulses could explain the observed 1/f 2 region in Figure 13, but this model alone is insufficient to incorporate the 1/f region. Let us now lift the requirement of equally decaying rates for bursts, and consider instead a system with a uniform distribution of decaying rates between two positive endpoints λ 1 and λ 2 , with λ 2 > λ 1 . The probability density function for a given λ in the interval is simply The PSD can still be calculated analytically by using this equation in combination with Equation (13), Equation (17) can be further simplified by considering the behavior of the arctan term between parenthesis in the domains 2πf ≪λ 1 , λ 1 ≪2πf ≪λ 2 and 2πf ≫λ 2 . The simplification in the first and last domains, respectively for very small and very large frequencies, can be done by linearizing the equation for small f and small 1/f . As for the domain with intermediate frequencies, the term between parenthesis is nearly constant and equal to π/2. The resulting behavior of S(f ) is We therefore conclude that the introduction of an uniform distribution of decaying rates for the pressure pulses is successful in predicting the 1/f to 1/f 2 transition in the power spectral density. This model also predicts the existence of a white noise region for very small frequencies (2πf ≪ λ 1 ) which is not present in our experiments. This is possibly because the value of λ 1 (that sets the upper limit in the frequency of this white noise domain) is too small to be noticeable in our setup. In the experiments, the smallest frequencies for which the 1/f domain holds are FIGURE 14 | Comparison between theory (thin blue line) and experiments (green crosses) for the CIP PSD. The theoretical prediction is given by Equation (19), in which we have used f t = 1.5 · 10 −2 Hz (marked as a dashed red line) and A 2 r = 1.35 Pa 2 /s. Guide-to-eye lines with exponents corresponding to the 1/f (pink noise) and 1/f 2 (brown noise) regions are also shown (thick black lines). The fit indicates the correctness of the asymptotic forms presented in Equation (20) and, more importantly, that this theory correctly captures the scaling features of the experimental PSD. Adapted from Moura et al. [23] under the Creative Commons CCBY license.
on the order of f = 7.0 10 −4 Hz, see Figure 13, which then means λ 1 < 2 π 7.0 10 −4 and for the sake of simplicity we can consider λ 1 = 0. Additionally, since λ 2 denotes the transition point between the 1/f to 1/f 2 regimes in Equation (18), we can write it as λ 2 = 2 π f t , where f t is the transition frequency. By plugging these assumptions into Equation (17), we are left with the simplified form which presents the asymptotic behavior From Figure 13 we can estimate the value of the transition frequency to be around f t = 1.5 · 10 −2 Hz for our experiment. By using this value, we can plot the analytical PSD Equation (19), to test the validity of the simplifications leading to it and, more importantly, to check how this analytical solution compares with our measurements. Figure 14 shows the resulting plot, together with our experimental measurements for the PSD and guide-toeye lines corresponding to the 1/f (pink noise) and 1/f 2 (brown noise) scaling regimes. The dashed red vertical line denotes the transition frequency f t . In producing this figure, we have used the product A 2 r as a tuning parameter for Equation (19), choosing the value A 2 r = 1.35 Pa 2 /s. We can see that 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 successfully, it also fits the data well for the crossover region between those domains. We therefore conclude that the modeling of the pressure bursts as a succession of decaying exponential pulses, with a uniform distribution of decaying rates successfully describes the scaling properties of the power spectral density for the CIP boundary condition. We remind that this modeling was motivated not only by the measured 1/f to 1/f 2 transition in the PSD, but also by the visual signature of bursts as decaying exponential pulses and their intermittent occurrence times, see Figure 7B. An additional explanation for the exponential behavior of the pressure bursts is given in Appendix 2 in Supplementary Material and a comment on the link between the pressure measured by the sensor and the actual capillary pressure across the liquid-air interface is found in Appendix 3 in Supplementary Material. Additional refinements of the model are still possible, even without losing the benefit of analytical tractability. One could, for instance, think of a model in which not only the pulses decaying rates λ, but also their amplitudes A follow a certain distribution. The calculations can be handled out in a similar manner, under the assumption that the distributions of λ and A are uncorrelated, the only difference being the fact that the constant value A 2 in Equation (17) would be exchanged by the expected value of this quantity. The resulting equation for a model with A being distributed according to a given probability density function ϑ(A) between two positive endpoint values A 1 and A 2 , with A 2 > A 1 and having λ again following a uniform distribution is in which the superscripts in S Ag denotes a general distribution in the amplitudes A, the distribution of the decaying rates λ remaining uniform. Notice that, this change will not interfere with the frequency dependency of S. Therefore, the power spectrum of a model with any distribution of amplitudes A and a uniform distribution of decaying rates λ still presents the same scaling behaviors as the one shown in Equation (18), in particular, it could also be used to explain the 1/f to 1/f 2 transition (pink to brown noise transition) observed in our experiments. Indeed, from Figure 7B, one would expect that the pulse amplitudes should follow some sort of distribution instead of assuming a single constant value. Nevertheless, considering that we are currently not measuring the product A 2 r (or its expected value in Equation 21) but merely using it as a fitting parameter in the theory, we have chosen to use the simplified version as in Equation (20) since this is sufficient to capture the relevant scaling features. Finally, as a last step of generalization in the model, one could consider the scenario in which the condition of a uniform distribution for the decaying rates λ is also lifted. By considering a power-law distribution of the form between the positive endpoints λ 1 and λ 2 , with λ 2 > λ 1 . The subscript in υ p denotes a power-law distribution and k is a normalization constant. The case δ = 0 corresponds to the one in Equation (14). Using this distribution, together with Equations (15) and (21) we obtain the general form for the power spectrum of a model with an arbitrary distribution ϑ(A) of amplitudes and a power-law distribution υ p (λ) of decaying times. This model is, surprisingly, still integrable (similar calculations can be found in Butz [93]). The integration of Equation (15) results in where, and 2 F 1 (a, b; c; z) is the Gauss hypergeometric function [90]. The subscript and superscript in S Ag λp denote a general distribution in the amplitudes A and a power-law distribution for the decaying rates λ. The fact that this model admits an analytical solution even for these conditions is surely worth noticing, nevertheless, the complicated form of Equations (23) and (24) does not yield much insight into the actual frequency dependency of S λp Ag . Since this dependency is entirely encoded in the I(f ) term in Equation (23), one can perform a linearization procedure to analyze the behavior for very small and very large frequencies, in a similar manner to what was done earlier to produce Equation (18). The resulting behavior for very small and very large frequencies is similar to what was obtained in the case of a uniform distribution of decaying rates, i.e., I(f ) ∝ 1/f 0 if 2πf ≪ λ 1 and I(f ) ∝ 1/f 2 if 2πf ≫ λ 2 . Therefore, the limiting cases exhibit the same white and brown noise profile. In the intermediate frequency range, i.e., λ 1 ≪ 2πf ≪ λ 2 , one can perform the change in variables u = λ/ 2πf to rewrite I(f ) as The resulting integral converges and we conclude that in the intermediate frequency range λ 1 ≪ 2πf ≪ λ 2 , I(f ) scales as I(f ) ∝ 1/f 1+δ . By plugging this result into Equation (23), we see that S Ag λp (f ) follows the same scaling, i.e.,

S
Ag By changing the exponent δ characterizing the distribution of decaying rates υ p (λ) this model can generate a whole family of 1/f 1+δ scaling regimes. Notice that the essential assumption in the model that is necessary to account for the 1/f region is the uniform distribution of decaying rates λ. Measuring this distribution directly from the pressure signal is not viable, given that recovering the underlying exponential bursts becomes very troublesome, particularly in the region in which several bursts appear in superposition (see inset of Figure 7). In Moura et al. [23] we show evidence for the claim that λ follows a nearly uniform distribution by considering the distribution of burst's durations (as defined in section 3.2).

CONCLUSIONS
We have performed an experimental investigation of two-phase flow phenomena in porous media, with particular focus on the question of how different boundary conditions influence the morphology of the flow and the underlying burst dynamics. We have observed that the invasion patterns are very similar for both boundary conditions, as long as the experiments are performed slowly (such as to guarantee that viscous forces are not relevant to the dynamics, which is dominated by capillary forces). Nevertheless, when it comes to dynamical features, some clear differences are observed. For example, the measured pressure signal presents pulses with a relaxation phase almost linear for the CWR boundary condition, whereas in the CIP case an exponential relaxation behavior is observed. We have analyzed the CWR pressure signal in terms of timedirected avalanches and showed how its scaling properties are in line with other works [85] in which those avalanches were investigated for systems with and without a correlated pore structure. Further on, we have shown that the differences between the CWR and CIP pressure signals are translated into their power spectral density, which, in the particular case of the CIP boundary condition, presents a 1/f scaling regime for low frequencies (pink noise), which then changes to 1/f 2 for intermediate frequencies (brown noise). This interesting scaling behavior is further described by a fully integrable mathematical framework in which the signal is modeled as a superposition of random exponential pulses. We have derived an expression for the PSD that was successful in describing our experiment (Equation 19) and showed how the model could be extended to a much more general class of systems, with less restricting assumptions on the distribution of pulse amplitudes and decay rates (Equations 23-26). By providing a rather general set of analytical expressions for the PSD, we hope this work will find applicability in a class of problems much wider than the specific case treated here.
Although all experiments performed employed quasi-2D samples, we believe some of our results could extend to 3D systems. In particular, the scaling features of the pressure signal power spectrum, should depend only on the intermittent dynamics of the pressure pulses, their exponential relaxation behavior and the distribution of decaying rates λ. We don't expect those properties to be a function of the dimensionality of the system, therefore, as long as the dynamics is driven in such a way that the system is in a quasi-equilibrium state, a CIP system in 3D should also present a similar 1/f to 1/f 2 transition in the PSD. Nevertheless, devising an experiment to perform this test in 3D is a challenging task, since the feedback mechanism used to guarantee the quasi-equilibrium state is based on the imaging of the flow and such imaging is harder to perform in 3D. One possible way to overcome this issue is to adapt the feedback mechanism to use the weight of the displaced fluid as a ruling parameter, thus lifting the need for full visualization of the flow.

DATA AVAILABILITY STATEMENT
The datasets related to this study can be obtained from the corresponding author upon request.

AUTHOR CONTRIBUTIONS
The experiments were performed and analyzed by MM, who was also responsible for writing the article. The experimental work was developed in the labs of KM, who also contributed to the definition of the research questions and the development of the experimental techniques employed. MM, KM, EF, and RT were involved in discussions that led to the theoretical analysis presented. All authors contributed actively in all stages of the work.

FUNDING
This work was funded by the Research Council of Norway through its Centre of Excellence funding scheme with project number 262644.

ACKNOWLEDGMENTS
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, the LIA France-Norway D-FFRACT, and the EU Marie Curie ITN FLOWTRANS network. We would also like to thank Mihailo Jankov, James Campbell, and Kristian Olsen for many fruitful discussions. Finally, we would like to express our gratitude for the very constructive feedback obtained from the reviewers in the production of this work.