Fredrik K Eriksen
email: f.k.eriksen@fys.uio.no

Renaud Toussaint

Antoine L Turquet

Knut J Måløy

Eirik G Flekkøy

Pneumatic fractures in confined granular media

Keywords: numbers: 83.60.Wc, 81.05.Rm, 47.20.Ma

We perform experiments where air is injected at a constant overpressure Pin, ranging from 5 to 250 kPa, into a dry granular medium confined within a horizontal linear Hele-Shaw cell. The setup allows us to explore compacted configurations by preventing decompaction at the outer boundary, i.e. the cell outlet has a semi-permeable filter such that beads are stopped while air can pass. We study the emerging patterns and dynamic growth of channels in the granular media due to fluid flow, by analyzing images captured with a high speed camera (1000 images/s). We identify 4 qualitatively different flow regimes, depending on the imposed overpressure, ranging from no channel formation for Pin below 10 kPa, to large thick channels formed by erosion and fingers merging for high Pin around 200 kPa. The flow regimes where channels form are characterized by typical finger thickness, final depth into the medium and growth dynamics. The shape of the finger tips during growth is studied by looking at the finger width w as function of distance d from the tip. The tip profile is found to follow w(d) ∝ d β , where β = 0.68 is a typical value for all experiments, also over time. This indicates a singularity in the curvature d 2 d/dw 2 ∼ κ ∼ d 1-2β , but not of the slope dw/dd ∼ d β-1 , i.e. more rounded tips rather than pointy cusps, as they would be for the case β > 1. For increasing Pin, the channels generally grow faster and deeper into the medium. We show that the channel length along the flow direction has a linear growth with time initially, followed by a power law decay of growth velocity with time as the channel approaches its final length. A closer look reveal that the initial growth velocity v0 is found to scale with injection pressure as v0 ∝ P 3 2 in , while at a critical time tc there is a cross-over to the behavior v(t) ∝ t -α where α is close to 2.5 for all experiments. Finally, we explore the fractal dimension of the fully developed patterns. For example for patterns resulting from intermediate Pin around 100-150 kPa, we find that the box-counting dimensions lie within the range DB ∈ [1.53, 1.62], similar to viscous fingering fractals in porous media.

I. INTRODUCTION

Several processes in engineering, industry and earth sciences involve pneumatic (gas) or hydraulic (liquid) fracturing of the soil, which occurs when fluids in the ground are driven to high enough pressures to deform, fracture and generate porosity in the surrounding soil or rock. For example in environmental engineering, pneumatic or hydraulic fracturing is done to enhance the removal of hazardous contaminants in the vadose zone (soil remediation) [START_REF] Suthersan | Remediation Engineering: Design Concepts[END_REF][START_REF] Goodroad | Potential of pneumatic fracturing to enhance in situ remediation technologies[END_REF], for soil stabilization injection to ensure a solid foundation for structures [START_REF] Kazemian | A review of stabilization of soft soils by injection of chemical grouting[END_REF], or in packer tests for project planning, risk assessment and safe construction of dams and tunnels [START_REF] Moayed | In-situ stress measurements by hydraulic fracturing method at gotvand dam site, iran[END_REF]. In industry, hydraulic fracturing is done to enhance oil and gas recovery [START_REF] Montgomery | Hydraulic fracturing: history of an enduring technology[END_REF][START_REF]Hydraulic fracturing whitepaper[END_REF][START_REF] Naik | Tight gas reservoirs -an unconventional natural energy source for the future[END_REF], CO 2 sequestration [START_REF] Johnson | Hydraulic fracturing mechanisms in carbon sequestration applications[END_REF], water well-and geothermal energy production [START_REF] Williamson | Hydraulic Fracturing to Improve the Yield of Bores in Fractured Rock[END_REF][START_REF] Rummel | The falkenberg geothermal frac-project: concepts and experimental re-sults[END_REF][START_REF] Clark | Energy company plans to frack volcano[END_REF]. Related natural processes such as subsurface sediment mobilization are studied in earth sciences, where sand injectites, mud diapirs and mud volcanoes are formed due to pore-fluid overpressure [START_REF] Hurst | Fluidization structures produced by upward injection of sand through a sealing lithology[END_REF][START_REF] Løseth | Gas and fluid injection triggering shallow mud mobilization in the hordaland group, north sea[END_REF][START_REF] Van Rensbergen | Near-surface sediment mobilization and methane venting in relation to hydrate destabilization in southern lake baikal, siberia[END_REF][START_REF] Talukder | Pliocene to recent mud diapirism and related mudvolcanoes in the alboran sea(western mediterrenean)[END_REF][START_REF] Pralle | Experimental evidence on the role of gas in sediment liquefaction and mudvolcanism[END_REF][START_REF] Deville | The origin and processes of mudvolcanism: new insights from trinidad[END_REF]. For example, the Lusi mud volcano in Indonesia is the biggest and most damaging mud volcano in the world [START_REF] Nuwer | Indonesias 'mud volcano' and nine years of debate about its muck[END_REF], having displaced 40 000 people from their homes, and has been active since May 2006. There is an ongoing debate about how it was triggered, i.e. whether it formed naturally by an earthquake or geothermal process [START_REF] Mazzini | Triggering and dynamic evolution of the lusi mud volcano, indonesia[END_REF][START_REF] Mazzini | Strike-slip faulting as a trigger mechanism for overpressure release by piercement structures. implications for the lusi mud volcano, indonesia[END_REF][START_REF] Mazzini | A new hydrothermal scenario for the 2006 lusi eruption, indonesia. insights from gas geochemistry[END_REF][START_REF] Lupi | Lusi mud eruption triggered by geometric focusing of seismic waves[END_REF], or that it is a man-made consequence of a nearby drilling operation by a company probing for natural gas [START_REF] Tingay | Initiation of the lusi mudflow disaster[END_REF].

Fluid injections into granular media has been extensively studied in laboratory experiments and simulations, where a common method to simplify the problem is to confine the experiment within a quasi-2-dimensional geometry, i.e. a Hele-Shaw cell. In [START_REF] Johnsen | Pattern formation during air injection into granular materials confined in a circular hele-shaw cell[END_REF][START_REF] Johnsen | Coupled air/granular flow in a linear hele-shaw cell[END_REF], the decompaction, fluidization regimes, and coupling between air and granular flow was studied in dry granular media in open circular and rectangular cells during air injection at different overpressures. Similar behavior was seen for liquid saturated granular material injected with the same liquid [START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas[END_REF], so it is reasonable to assume that studies of pneumatic fracturing also have qualitative relevance to hydraulic fracturing. The two cases are similar in the sense that no interface tension and viscosity difference exist between the invading and interstitial fluids. In systems where gravity is felt significantly, as in inclined or vertical systems, the density of the fluid is a notable difference, and buoyancy forces on the grains are significant for liquids, but usually negligible for gases. In horizontal layers and thin compacted systems, where gravity is perpendicular to the system and plays no major role on the horizontal motion of such quasi two-dimensional systems, the main difference in an experiment where a liquid is used instead of a gas is that the compressibility is lower such that the overpressure imposed at the boundary diffuses much faster [START_REF] Niebling | Sedimentation instabilities: impact of the fluid compressibility and viscosity[END_REF]. However, near the invading channels, the pressure gradients are found to be qualitatively similar in both cases [START_REF] Eriksen | Pressure evolution and deformation of confined granular media during pneumatic fracturing[END_REF], leading to similar growth mechanisms for the channels [START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas[END_REF]. Further, the compressibility felt in water on the scale of a ∼ 100 m reservoir is similar to the compressibility felt in air on the scale of a ∼ 1 m experiment, and therefore the effect of pressure diffusion is included when air is used in ∼ 1 m model reservoirs. Another difference is that for a given pore pressure gradient and grain configuration (permeability), the fluid velocity between the grains is inversely proportional to the viscosity force. The stress thresholds allowing to overcome friction with the confining plates depends mainly on the solid fraction, and does not change in order of magnitude between cases where the interstitial fluid is a gas or a liquid [START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas[END_REF]. Consequently, for similar channel geometries and grain packings, the pressure difference between inlet and outlet that lead to grain motion are similar in the two systems. However, the grain mobility is observed to be essentially inversely proportional to the pore fluid viscosity, i.e. the observed dynamics are similar in terms of a reduced time corresponding to time multiplied by the pore fluid viscosity [START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas[END_REF]. The patterns observed are also close to each other [START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas[END_REF]. The patterns formed during fluid injection into a granular medium, and evolution of the fluid-solid interface, have been found to resemble Diffusion Limited Aggregation patterns (DLA) and viscous fingering [START_REF] Cheng | Towards the zero-surface-tension limit in granular fingering instability[END_REF], a fingering instability that occurs when a less viscous liquid is injected into a porous medium containing a more viscous liquid with which it cannot mix [START_REF] Saffman | The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid[END_REF]. As mentioned in [START_REF] Johnsen | Coupled air/granular flow in a linear hele-shaw cell[END_REF], the main difference between the viscous-and granular fingering instabilities is the absence of interfacial tension in the granular case. For example, the stabilizing forces in viscous fingering are surface forces, while in granular fingering it is the build-up of friction between particles and against the confinement. However, both instabilities are driven by the pressure gradient across the defending medium, which is largest on the longest finger tips, making more advanced fingers grow on expense of the less advanced ones. A notable difference between air injection into a dry granular medium and a saturated one is that the overpressure ini-tially diffuses into the packing in the dry case, while it is already a steady-state Laplace field over the defending liquid in the saturated case. This difference arises from the compressibility of the defending phase.

Further, during air injection into liquid saturated granular media and suspensions, the characteristics of emerging patterns and behavior of the media depend on injection rate, and the competition between mobilized friction and surface forces [START_REF] Eriksen | Invasion patterns during two-phase flow in deformable porous media[END_REF][START_REF] Sandnes | Patterns and flow in frictional fluid dynamics[END_REF][START_REF] Holtzman | Capillary fracturing in granular media[END_REF][START_REF] Kong | Morphodynamics during air injection into water-saturated movable spherical granulates[END_REF][START_REF] Chevalier | Morphodynamics during air injection into a confined granular suspension[END_REF][START_REF] Eriksen | Bubbles breaking the wall: two-dimensional stress and stability analysis[END_REF][START_REF] Eriksen | Numerical approach to frictional fingers[END_REF][START_REF] Marks | Compaction of granular material inside confined geometries[END_REF][START_REF] Eriksen | Pattern formation of frictional fingers in a gravitational potential[END_REF][START_REF] Moura | Impact of sample geometry on the measurement of pressure-saturation curves: Experiments and simulations[END_REF][START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF][START_REF] Løvoll | Growth activity during fingering in a porous hele-shaw cell[END_REF][START_REF] Løvoll | Influence of viscous fingering on dynamic saturationpressure curves in porous media[END_REF][START_REF] Tallakstad | Steadystate two-phase flow in porous media: Statistics and transport properties[END_REF]. For example, one observes flow regimes such as two phase flow in rigid porous media [START_REF] Moura | Impact of sample geometry on the measurement of pressure-saturation curves: Experiments and simulations[END_REF][START_REF] Wilkinson | Invasion percolation: a new form of percolation theory[END_REF][START_REF] Løvoll | Growth activity during fingering in a porous hele-shaw cell[END_REF][START_REF] Løvoll | Influence of viscous fingering on dynamic saturationpressure curves in porous media[END_REF][START_REF] Tallakstad | Steadystate two-phase flow in porous media: Statistics and transport properties[END_REF], capillary fracturing, stick-slip bubbles and labyrinth patterns [START_REF] Eriksen | Invasion patterns during two-phase flow in deformable porous media[END_REF][START_REF] Sandnes | Patterns and flow in frictional fluid dynamics[END_REF][START_REF] Holtzman | Capillary fracturing in granular media[END_REF][START_REF] Kong | Morphodynamics during air injection into water-saturated movable spherical granulates[END_REF][START_REF] Chevalier | Morphodynamics during air injection into a confined granular suspension[END_REF][START_REF] Eriksen | Bubbles breaking the wall: two-dimensional stress and stability analysis[END_REF][START_REF] Eriksen | Numerical approach to frictional fingers[END_REF][START_REF] Marks | Compaction of granular material inside confined geometries[END_REF][START_REF] Eriksen | Pattern formation of frictional fingers in a gravitational potential[END_REF]. In the opposite case, during liquid injection into dry granular media [START_REF] Huang | Granular fingering in fluid injection into dense granular media in a hele-shaw cell[END_REF], for a given imposed flux, the flow behavior goes from stable invasion towards saturated granular fingers for increasing flow rate and viscosity of the invading fluid. These fingering patterns are thought to form due to the permeability contrast between the channels empty of grains and the granular medium, and are in this sense similar to viscous fingering, with a lower viscous pressure drop in empty channels than in the porous medium, as in classical viscous fingering systems where the viscosity contrast leads to a similar effect, i.e. where the invading fluid has a lower viscosity than the invaded one [START_REF] Saffman | The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid[END_REF][START_REF] Løvoll | Growth activity during fingering in a porous hele-shaw cell[END_REF][START_REF] Toussaint | Influence of pore-scale disorder on viscous fingering during drainage[END_REF]. An important difference with classical viscous fingering is the presence of solid stresses and solid friction in the granular material. This is a stabilizing phenomenon absent from systems with two Newtonian fluids: A threshold arising from friction with the system boundaries has to be overcome in viscous fluid drag to allow for grain motion. The same trend is shown in numerical studies for gas injection into granular media containing the same gas [START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF].

Granular fingering instabilities have also been studied in closed vertical cells, where gravity drives the flow as heavier beads fall down from a granular layer at the top of a lighter fluid layer [START_REF] Niebling | Mixing of a granular layer falling through a fluid[END_REF][START_REF] Vinningland | Granular rayleigh-taylor instability: Experiments and simulations[END_REF][START_REF] Vinningland | Experiments and simulations of a gravitational granular flow instability[END_REF][START_REF] Vinningland | Size invariance of the granular rayleigh-taylor instability[END_REF][START_REF] Vinningland | Family-vicsek scaling of detachment fronts in granular rayleightaylor instabilities during sedimentating granular/fluid flows[END_REF][START_REF] Völtz | Rayleigh-taylor instability in a sedimenting suspension[END_REF]. When the beads detach at the front, they form fingers of falling granular material surrounding finger-like bubbles of rising fluid. These fingers are found to coarsen over time until they reach a typical wavelength depending on the interstitial fluid and bead size.

Typically, in all processes involving fluid injection into granular media, there are flow regimes where the medium has either solid-like behavior or fluid-like behavior. This is one of the special properties of granular materials, which also show gas-like behavior in some cases [START_REF] Jaeger | Granular solids, liquids, and gases[END_REF].

In this paper, we present an experimental study on flow regimes and pattern formation during air injection into confined granular media, while the surrounding deformation and pore-pressure evolution is studied in a related article [START_REF] Eriksen | Pressure evolution and deformation of confined granular media during pneumatic fracturing[END_REF]. More specifically, we inject air at constant overpressure into a dry granular medium inside a Hele-Shaw cell, where air escapes at the outlet while beads cannot. The motivation of this setup is to study the granular Saffman-Taylor instability in compacting granular media, and the coupling between compaction and flow. A similar, but smaller system has been studied in numerical simulations by Niebling et al. [START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF][START_REF] Niebling | Numerical studies of aerofractures in porous media[END_REF]. During air injection at different overpressures they found that fractures grow faster, longer, as well as coarsen with increasing injection pressure. There, the fracture propagation velocity is roughly constant initially and scales with the square root of the injection pressure. Further, by varying the interstitial fluid viscosity, two flow regimes were identified; one with finely dispersing bubbles and large scale collective motion of particles, the other one with build-up of a compaction front and fracturing. These flow regimes depend, respectively, on whether the particles are primarily accelerated by the imposed pressure gradient in the fluid, or interactions through particle contacts. This in turn depends on the diffusivity of the interstitial fluid pressure in the granular medium. We analyze our findings in light of these results, to investigate what is similar and what is different in our roughly 10 times larger cell with the same cell gap. As opposed to similar experiments, with open outer boundary conditions, after the flow compacts the medium there is no decompaction. We thus expect the material behavior (at high enough overpressure to displace beads) to have a transition from fluid-like to solid-like during experiments, and that eventual invasion patterns will initially resemble viscous fingering in the fluid-like regime, crossing over to stick-slip fracture propagation as the medium becomes more solid-like, until it reaches a final structure as the compacted medium has reached a completely solid-like behavior. What is less obvious, is how the flow patterns in this system change with the injection pressure. By varying the imposed overpressure, we identify and describe the different flow regimes. Due to the confined nature of the experiment, it is thought to be a laboratory analog to pneumatic and hydraulic fracturing of tight rock reservoirs where the free boundary at the surface is very distant from the injection zone. Therefore, new insight into this problem may have industrial applications in addition to increase the understanding of flow and transformations in porous media.

It is also worth to mention a closely related project [START_REF] Turkaya | Bridging aero-fracture evolution with the characteristics of the acoustic emissions in a porous medium[END_REF], where acoustic emissions recorded during the experiments are analyzed. There, it is shown that different stages of the invasion process can be identified acoustically in terms of characteristic frequencies and distinct microseismic events.

II. METHODS

A. Experimental setup

The experimental setup is a linear Hele-Shaw cell, partially filled with Ugelstad spheres [START_REF] Ugelstad | Swelling of oligomer-polymer particles. new methods of preparation[END_REF], i.e. dry, nonexpanded polystyrene beads with a diameter of 80 µm ± 1 %. The cell is made out of two rectangular glass plates (80×40×1 cm in length, width and thickness respectively) clamped together on top of each other with an aluminum spacer controlled separation of 1 mm. A cell volume (76×32×0.1 cm) is formed between the plates by an impermeable sealing tape as shown in figure 1, with one of the short sides left open (outlet). Next, beads are filled into the cell by pouring them through the open side until the packing occupy about 90 % of the cell volume, followed by closing the open side with a semi-permeable filter (a 50 µm steel mesh) to keep beads inside the cell while allowing air to escape. The cell is then flipped vertically to place the granular medium against the semipermeable outlet by using gravity, resulting in a volume packing fraction of approximately ρ s = 0.44 ± 0.04, assumed to be more or less uniform across the medium (in a quasi-2D system of ∼ 10 beads thickness, the confining plates cause a frustrated system with lower packing fraction [START_REF] Nelson | Order, frustration, and two-dimensional glass[END_REF], as observed previously in similar systems where ρ s = 0.44 [START_REF] Johnsen | Pattern formation during air injection into granular materials confined in a circular hele-shaw cell[END_REF] and ρ s = 0.42 [START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF]). This leaves a volume empty of beads on the sealed side of the cell, opposite to the semi-permeable outlet, with a linear air-solid interface. An inlet hole on the sealed side of the cell is connected to a pressurized air tank which lets us inject air at a constant and maintained overpressure, P in = P abs,in -P 0 (absolute pressure -atmospheric pressure), ranging from 5 to 250 kPa. This will force air to move through the granular medium, towards the semi-permeable outlet, where P abs,out = P 0 = 100 kPa, or in terms of overpressure above the atmospheric one, P out = 0.

During experiments, the prepared cell is positioned horizontally. A selected overpressure is set at the pressure tank outlet and verified by a Honeywell pressure sensor with an accuracy of ±4 kPa. The tubing between the pressure tank and the cell inlet is equipped with an electronic valve such that the air injection is started with a digital trigger signal. Positioned above, with a topdown view of the cell, a Photron SA5 high speed camera is started with the same trigger signal, recording the air invasion at a framerate of 1000 images/s and a resolution of 1024×1024 pixels (1 pixel ≈ 0.7 mm in the cell). Light from a 400 W Dedolight studio lamp provides uniform and flicker-free illumination onto the white beads of the medium. A small fraction of the beads (< 10 %) are dyed black with ink to create tracer particles that are used for tracking frame-to-frame deformation in the granular medium, which is studied in another article [START_REF] Eriksen | Pressure evolution and deformation of confined granular media during pneumatic fracturing[END_REF] (note that in figures 2, 3 and 6, there are some darker areas creating patterns in the initial packing. This is due to a higher concentration of dyed beads, which is related to the initial preparation process and the packing method. It does not seem to have any influence on the experiments). The experiments are run for 10 s, but typically the fracturing and/or compaction of the granular medium takes less than 5 s.

B. Image processing

In analysis of the images from the high speed camera, we investigate flow regimes, characteristics and evolution of the fingering patterns formed. We perform image processing with Matlab to obtain the information contained in the images. Quantitative analysis of the invasion patterns is done by converting the grayscale raw data into binary images, i.e. images with either black or white pixels, where the white pixels (value=1) represent the pattern and the black pixels (value=0) represent the background, or the bead-filled region. We obtain such binary image sequences from the raw data by thresholding each frame with the initial image, such that the pixels having a value less than 30% of the corresponding initial value becomes white and the rest remain black, as shown in figure 2 (c). In addition, once a pixel is invaded (and made white) it will remain white for the rest of the experiment to correct for mobilized beads flowing inside channels, which occurs due to erosion and fingers merging in some experiments. From the binary images we extract information such as invasion depth over time, average finger thickness as function of depth, and fractal dimensions.

III. RESULTS

Depending on injection pressure (and friction related to packing fraction) we observe the granular medium to exhibit either a solid-like behavior, or a transition from fluid-like to solid-like behavior during air injection. In the solid-like regime, there is no apparent deformation (bead displacements), and the air is reaching the cell outlet by seeping through the network of pores between beads. On the other hand, in the fluid-like regime we observe significant deformation, where beads are displaced by an amount corresponding to several bead sizes, by momentum exchange between air molecules and beads. At sufficient overpressure the granular medium has a behavior much like a viscous liquid being invaded by air, as the flow opens up channels empty of beads, a Saffman-Taylor like instability [START_REF] Saffman | The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid[END_REF] resulting in dendritic invasion patterns. A stabilizing mechanism of this flow instability is the build-up of friction and stress chains between beads during compaction, so due to the boundary conditions in our experiments we always end up with a solid-like medium towards the end of an experiment. Since the channels empty of beads display patterns with finger-like features, and that the air-solid interface expands similarly to a fingering instability, we refer to the branches of the channels as fingers when we characterize the patterns. At the same time, the fingers have a long and slender morphology which is fracture-like, and during growth frictional "bonds" between the compacting grains and the confining plates are broken in order to open the channels. Therefore we consider the channel growth as fracturing, although there are no cohesive bonds broken as in the traditional sense of fracturing.

A. Flow regimes

To explore the various flow regimes we have systematically varied the air injection pressure within the range of 5 kPa to 250 kPa, and we have identified 4 qualitatively different flow regimes occurring in this range. Within certain flow regimes, the pattern characteristics gradually change with increasing injection pressure, so we have grouped the experimental results into six categories defined by final pattern characteristics and how the channels expand. The categories are chosen in order to obtain average results from similar experiments, and are defined in the list below. Typical final structures are shown for each category in figure 3, and the associated range of injection pressure is listed in table I. Videos of channels formed in categories 3a,b and 4a,b are found in the supplementary material [60-63].

• Category 1 -No deformation:

The porous medium appears rigid (solid-like), with no significant bead displacement or invasion pattern formed. This is flow regime #1, where the pore pressure gradient does not overcome the threshold for bead displacement which arise from granular stress, friction between the beads and bead-plate friction.

• Category 2 -Initial compaction without channels:

The medium is initially compacted during a collective displacement of beads towards the outlet side, without the formation of invasion channels. This is flow regime #2, where the pore pressure gradient is initially high enough to overcome the granular stress threshold arising from the friction with the plates. Then, as the pore pressure diffuses into the medium, the pore pressure gradient drops below the bead displacement threshold and the deformation stops before channels are formed.

During the short time beads are displaced, they are primarily moved by the pore pressure gradient, i.e. stress from the fluid injection transmitted through solid contacts (compacting zone) does not significantly extend beyond the influence of the overpres-sure (skin depth), leading to collective bead displacement rather than the formation of fractures [START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF].

• Categories 3a and 3b -Thin branched channels are formed: After the initial compaction, a dendritic invasion pattern is formed with an average finger thickness less than 2 cm. This is flow regime #3, where the pressure gradient is high enough to deform the medium, also after the compacting zone reaches beyond the influence of the pore pressure gradient. This leads to random perturbations in the airsolid interface and fracture channels begin to form [START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF]. At later times, the channel growth stops as the medium is compacted, i.e. the solid stress and friction increase such that the bead displacement threshold eventually becomes too large for further deformation. In this flow regime, the channels become longer with increasing injection pressure, so we define arbitrarily (in order to perform statistical analyses on groups of experiments) category 3a as when the final length of the channel is less than 50% of L = 70 cm, the initial length of the granular medium, and category 3b as when the final length of the channel is more than 50% of L.

• Categories 4a and 4b -Branched channels are formed and eroded: After the initial compaction, a long dendritic invasion pattern is formed (as in category 3b) with a main channel thicker than the branch fingers, due to erosion inside the channel. This is flow regime #4, where the high pressure gradients at the most advanced channel tips cause a large enough fluid flux inside the channels to re-mobilize beads along with the flow. The width of the main channel increases by erosion (i.e. motion of beads nonperpendicular to the channel interface, but with a major tangential component along it). For higher pressures, side channels and branches may also merge into the main channel. In this case, the high flow velocity in the main channel might reduce the fluid pressure slightly such that adjacent channels grow into it. In this flow regime, the degree of erosion increases gradually with increasing injection pressure, so we define two categories. Category 4a is defined as when there is only erosion, where the average finger thickness of the final structure is typically between 2 cm and 4 cm. The erosion is not persistent, and rarely occurs more than 15 cm behind the most advanced tip. Category 4b is defined as when there is both erosion and channels merging, where the average finger thickness of the final structure is typically more than 4 cm. Erosion and channel merging is observed to occur even at distances of 30-40 cm behind the most advanced tip. Again, the erosion is not persistent and stops around the point when the channel structure reaches 90 % of its final length.

As one can expect due to the heterogeneous nature of the grain packings, and the variability of the grain assembly between experiments, a significant dispersion of the resulting patterns is observed for experiments at similar imposed pressures, and consequently, the pressure ranges for the categories overlap (see table I). The distribution of new growth dA behind the longest fingers at different intervals is shown in figure 4 for typical experiments in the categories 3a,b and 4a,b. The new growth dA refers to the change in channel area between two successive snapshots. The new growth in figure 4 is found by subtracting successive binary images A of the patterns pixel by pixel, i.e. dA = A(t + ∆t) -A(t), where ∆t is the time between snapshots where the channels grow by 10 % increments of their final lengths (this method of selecting ∆t is done in order to determine where the channels are growing at different stages, and thus to plot distributions of dA with comparable space integrals at different times of the experiment, since the growth rate decreases at later times). In the plots, the distribution of new growth over the distance d behind the most advanced finger is represented by a normalized value n dA (d), which is the mass of the new growth dA at the depth d, divided by the width of the cell such that n dA (d) is a dimensionless number ∈ [0, 1] indicating the fraction of the cell width at position d that contains new growth. For categories 3a,b and 4a the growth typically extends in a region of 10-15 cm behind the most advanced tip, whereas for category 4b the growth extends along the whole channel length. The category 4b channels have a rapid and large volume change compared to the other ones, with a high air flow velocity resulting in erosion along the entire channel. In addition, the rapid expansion could let the channel pressure remain somewhat lower than P in during growth, such that when it stops the pressure rises to P in . If it is so, pressure gradients could form at the solid-air interface far behind the most advanced finger explaining the small growth there. The displacement of the granular medium outside the channels is similar for all categories; Beads are pushed in front of the channels and compacted until the system is jammed. Beads are also pushed out from the side of the advancing channel tips, opening up for branches.

In the rest of the paper, we focus on the categories 3a,b and 4a,b, which are in flow regimes #3 and #4 where invasion channels are formed. For these flow regimes, we discuss the typical characteristics of the patterns formed and their growth dynamics. Cat. 3a The new growth is measured on intervals between successive snapshots where the channel length grows 10 % of the final channel length for each snapshot. The large initial new growth (0 % -10 %) in all plots corresponds to the initial bulge that forms before the channeling instability begins. The plots indicate the depth of the active growth zone, which goes from very narrow for category 3a (the growth is focused on the most advanced tips), to the whole channel length for category 4b.

Cat. 3b

LocalLgrowthL/LcellLwidth,Ln dAL 0L-L10LA 10L-L20LA 20L-L30LA 30L-L40LA 40L-L50LA 50L-L60LA 60L-L70LA 70L-L80LA 80L-L90LA 90L-L100LA 0L-L10LA 10L-L20LA 20L-L30LA 30L-L40LA 40L-L50LA 50L-L60LA 60L-L70LA 70L-L80LA 80L-L90LA 90L-L100LA Cat. 4a
0L-L10LA 10L-L20LA 20L-L30LA 30L-L40LA 40L-L50LA 50L-L60LA 60L-L70LA 70L-L80LA 80L-L90LA 90L-L100LA Cat. 4b 0L-L10LA 10L-L20LA 20L-L30LA 30L-L40LA 40L-L50LA 50L-L60LA 60L-L70LA 70L-L80LA 80L-L90LA 90L-L100LA DistanceLbehindLfingerLtip,LdL[cm]

B. Finger characteristics

In figure 5 we see the average finger thickness w of the final structures as function of depth x into the granular medium, measured from the initial boundary position. The average finger thickness of a pattern at a given depth into the granular medium is found by intersecting the binary image of the pattern with a line perpendicular to the average flow direction at that depth, and count the total number of white pixels intersected. In other words, the total width of the pattern is measured followed by dividing it by the number of connected pixel groups, i.e. the number of fingers. This is done for each experiment and the results are averaged within each category. We see that the patterns in categories 3a and 3b have roughly the same thickness, averaging around 1.4 cm and is more or less constant as function of depth. We do not observe erosion inside these channels, so a finger thickness around 1.4 cm seems typical when there is no erosion. The category 4a patterns show initially thicker fingers along the first half of the structures, here with an average thickness of 3.3 cm, which then decreases down towards the typical thickness seen for categories 3a and 3b. The larger initial finger thickness is due to erosion inside the main channel, meaning that beads near the channel walls are re-mobilized in the flow direction after the initial growth, making the channel grow perpendicular to it. As average thicknesses here suggest, this effect typically increase the channel width by an amount of around 2 cm. The category 4b patterns typically have thick fingers throughout their entire length up to the tips, with an average of 5.3 cm. This thickness is achieved both by erosion inside the main channel, and with fingers merging together. After fingers merge and trap clusters of beads, these beads are mobilized in the flow direction. In general we see an increase in finger thickness and length for increasing pressures. Some examples of both erosion and finger merging are shown in figure 6 for an experiment with P in = 200 kPa.

The plot in figure 7 shows the average number of fingers N f as function of depth x into the medium (found as the number of connected pixel groups on vertical crosssections, as explained for figure 5) for the final structures in each category. For patterns in the categories 3a and 3b, the number of fingers increase at lower depths than for patterns in the categories 4a and 4b, where there is a more pronounced initial region with fewer fingers, crossing over to more fingers during the last half of their length. This could indicate that channels with few branches form when the friction is low compared to the driving force ∝ ∇P , with a cross-over to the patterns branching out with more fingers as the friction is becoming comparable to the driving force due to compaction. The sudden drop in N f at the deepest parts of the patterns is due to the finite size of the structures, where only a few of the longest fingers reach.

To get a description of the shape of finger tips, we look at how the finger thickness w grow as function of distance d from the finger tip, which is found to follow a power law within the length of the tip region [START_REF] Cheng | Towards the zero-surface-tension limit in granular fingering instability[END_REF] (where the finger thickness increase with d before it saturates): (

) 1
When w is plotted as a function of d in a log-log plot, the slope β for w(d) within the tip region reveal information about the shape of the tip; it has a more rounded shape for β < 1, while it has a pointy cusp-shape for β > 1. We have averaged w(d) over time for each experiment, i.e. obtaining the average finger tip shape during flow in each experiment. Figure 8 (a) shows a log-log plot of the finger thickness as function of distance from the finger tip, averaged per category. The collapsed average slopes from all experiments seem to fall along the same line with β = 0.68, indicating more rounded finger tips with the same shape in all experiments. However, a purely round profile would give β = 0.5, so the fact that we measure a bigger value means that there still is a singularity in curvature towards the tip -even if it is not a spike shape. The singularity can be shown as follows, where

w ∼ d β ⇒ d ∼ w 1 β ⇒ dd dw ∼ w 1 β -1 , (2)
gives the slope of the front. From this we find the expression for the curvature as,

κ ∼ d 2 d dw 2 ∼ w 1 β -2 (or d 1-2β ). (3)
If β > 0.5, 1 β -2 < 0 and κ diverges at the tip where w (and d) → 0. At the same time, if β < 1, the slope does not diverge at d = 0, and yields a more rounded profile.

We have not observed any cusps as reported in [START_REF] Cheng | Towards the zero-surface-tension limit in granular fingering instability[END_REF] with β = 1.43 ± 0.2, consistent with theory on fluid fingering in the zero-surface tension limit. The evolution of β over time is plotted in figure 8

C. Fractal analysis

By looking at the invasion patterns, we notice that they have self-similar features, e.g. a smaller branch resemble the whole larger pattern. To characterize this feature we analyze fractal dimensions of the final patterns in three different ways. The box-counting dimension D B is found for each final structure by covering the binary image with boxes, i.e. dividing the image into equal squares of sides s, and count the number N of squares that contain a white pixel as function of box size s. For a fractal, the relationship between the number of boxes covering the pattern and their size follow a power law

N (s) ∝ s -D B , (4)
such that the box-counting dimension D B is found as the negative value of the slope of N (s) in a log-log plot [START_REF] Eriksen | Invasion patterns during two-phase flow in deformable porous media[END_REF][START_REF] Toussaint | Influence of pore-scale disorder on viscous fingering during drainage[END_REF][START_REF] Barabási | Fractal Concepts in Surface Growth[END_REF][START_REF] Feder | Fractals[END_REF]. By obtaining box-counting data over a range of sizes s, we estimate D B from the slope of linear fits between an upper cutoff s = 32 cm (cell width) and a lower cutoff s = 1 cm (typical for thinner fingers).

In addition, we obtain the box-counting dimension D F of the front (air-solid interface). This is found in the same way as the box dimension, but with binary images where only pixels on the perimeter of the patterns are white. Finally, we estimate local fractal dimensions D L as function of depth along the structures. To get the local fractal dimension at a given depth, we intersect the structure at that depth with a vertical line (perpendicular to the flow direction) and do a 1-dimensional box counting along that line, i.e. divide the line into pieces of equal length l and count the number N of line segments containing white pixels as function of l. Again, for a fractal we have the power law

N (l) ∝ l -D L , (5)
and we find the local fractal dimensions at given depths from slopes of linear fits between l = 16 cm and l = 1 cm. To compare D L with D B we use one of Mandelbrot's rules of thumb [START_REF] Feder | Fractals[END_REF][START_REF] Mandelbrot | The fractal geometry of nature[END_REF]. It states that the codimension of an intersected set equals the sum of the codimensions of the individual intersecting sets, here given by

E 2 -D L = (E 2 -E 1 ) + (E 2 -D) ⇓ D = D L + 1, ( 6
)
where D is the fractal dimension of the pattern, E 1 = 1 is the dimension of the line intersecting it, and E 2 = 2 is the dimension of the image plane containing the sets.

The box-counting dimensions are averaged per flow category and plotted in figure 9. We see that the box dimension for the less developed category 3a patterns is D B = 1.41. For categories 3b and 4a, the box dimensions are D B = 1.54 and D B = 1.53 respectively. For category 4b patterns, we observe D B = 1.63. Patterns in category 3b and 4a have box dimensions within the range of fractal dimensions found for viscous fingers in saturated porous media, i.e. D = 1.53 -1.62 [START_REF] Eriksen | Invasion patterns during two-phase flow in deformable porous media[END_REF][START_REF] Toussaint | Influence of pore-scale disorder on viscous fingering during drainage[END_REF][START_REF] Måløy | Viscous fingering fractals in porous media[END_REF]. Figure 10 shows the results of box-counting the airsolid interface of the final structures, averaged per category. We see that categories 3a,b and 4a have interface dimensions D F = 1.40, 1.52, 1.48 respectively which are close to their corresponding box-counting dimensions, while category 4b patterns have D F = 1.41 which is about 0.2 less than the corresponding box dimension. This can be interpreted in terms of finger thickness, i.e. the front becomes more similar to the structure for decreasing finger thickness.

Figure 11 shows the local fractal dimensions D L as function of depth into the granular medium. For all flow categories, the curves show a more or less constant behavior over most of the depths (the initial high dimension is due to a larger initial area empty of beads, and the final lower dimension is due to finite size effects of the patterns), which indicates that the patterns have consistent local fractal dimensions over a range of x. The categories 3a,b and 4a patterns have typical local dimensions of 1.53, 1.58, 1.60, respectively. This is within the range of dimensions seen for viscous fingers in saturated porous media (D = 1.53 -1.62), while for the category 4b patterns the local dimension is higher with a typical value of 1.76, which is closer to the dimensions seen for DLA clusters and viscous fingers in an empty saturated Hele-Shaw cell (D = 1.71).

D. Growth dynamics

In figure 12 (a), the finger tip position averaged per category is plotted as function of time, which is found by recording the maximum depth of the patterns in the binary images at each snapshot. We see that for all experiments, the channel length grows linearly with time initially before it decelerates towards a final length, and that both the growth velocity and final depth of the invasion structures increase with increasing injection pressure. The initial constant growth rate seems to scale with the injection pressure as P in to show that they fall along the same line initially. In the simulations done by Niebling et al. [START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF][START_REF] Niebling | Numerical studies of aerofractures in porous media[END_REF] of a similar, but smaller system they too observe an initial linear growth velocity crossing over to decay over time, however we here observe the initial growth velocity to scale with P in as v ∝ P 3/2 in instead of v ∝ P 1/2 in as found in these simulations. In figure 13 (a), the log-log plot of the finger tip velocity v(t)/(P 3/2 in ) as function of time indicates that when the finger tip velocity begins to decrease, it follows a power law equation with time, v(t) ∝ t -α , where α = 2.5 ± 0.2. in , where Pin = 64, 110, 134 and 213 kPa, for category 3a,b and 4a,b respectively. The initial constant plateau corresponds to the period with linear growth over time, which then crosses-over to a declining slope close to -2.5 (dashed line) for all the categories. The slopes show that the decay in growth velocity follows a power law equation with time. (b) The curves of normalized finger tip velocity v(t * )/v0 are collapsed along the time axis by scaling the time as t * = t • (v0/x f ) = R • (t/tc), where x f is the final channel length. The knee point of the curves, R ≈ 0.6, is marked by the dashed vertical line.

Timeb[s]

E. Scaled growth model

As we have shown in figures 12 and 13, the growth of the channel length is linear with time initially until it crosses over to decay towards a final channel length. More specifically, the growth velocity is constant initially, before it crosses over to a power law decay with time. The constant initial velocity v 0 is found to scale with injection pressure P in as

v 0 = C • P 3 2 in , (7)
where C ≈ 10 -3 2 cm/(s•kPa

2 ) is found from the initial part of the collapsed curves v(t)/P 3 2 in along the ordinate axis in figure 13 (a). By inserting the value for C into equation [START_REF] Naik | Tight gas reservoirs -an unconventional natural energy source for the future[END_REF], we can instead collapse the curves along the ordinate axis with v(t)/v 0 resulting in dimensionless and normalized units for the velocity. In order to collapse the curves along the time axis, we define a critical time t = t c as the time when the growth velocity crosses over from being constant to follow a power law decay with time. At the critical time, the finger tip position can be described mathematically as

x c = x(t c ) = v 0 • t c , (8)
which gives the relation

t c = x c v 0 = x c C • P 3 2 in ( 9
)
for the critical time.

Since we can easily measure x f , the final length of the channel, we check if x c can be described as a typical fraction R of x f , such that x c = R • x f , which gives

t c = R • x f C • P 3 2 in = R • x f v 0 . (10)
By plotting v(t * )/v 0 as function of

t * = t • (v 0 /x f ) = R • (t/t c
) in a log-log plot, we do indeed find that the curves collapse along the time axis with a crossover point corresponding to R ≈ 0.6 (read from the knee point in figure 13 (b)). Since we now have the constants C and R, we can calculate v 0 and t c for individual experiments by inserting the respective P in and x f into equations [START_REF] Naik | Tight gas reservoirs -an unconventional natural energy source for the future[END_REF] and [START_REF] Rummel | The falkenberg geothermal frac-project: concepts and experimental re-sults[END_REF]. Figure 14 shows log-log plots of v = v(t )/v 0 as function of t = t/t c for individual experiments with injection pressures in the range of 50 -250 kPa, which follow the same dimensionless curve given by the function

f v (t ) = 1, if t ≤ 1 t -α , if t > 1, where α = 2.5. (11)
Similarly, we collapse the fingertip position x = x(t )/x c and plot it as function of t = t/t c for the same individual experiments, also shown in figure 14 (Note that the data is from single experiments, i.e. not the average values as in figure 13, explaining the increased amount of noise). Here, the collapsed data follow the curve given by the function

f x (t ) = t , if t ≤ 1 1 1-α t 1-α + α α-1 , if t > 1, (12)
which is obtained by integrating each part of equation [START_REF] Clark | Energy company plans to frack volcano[END_REF], requiring that f x (t > 1) → 1 when t → 1. Further, f x (∞) = α/(α -1) = x f /x c = 1/R, giving R = (α -1)/α = 0.6, which fits very well with our observation in figure 13 (b). By substituting α into equation [START_REF] Rummel | The falkenberg geothermal frac-project: concepts and experimental re-sults[END_REF] we get,

t c = α -1 α • x f C • P 3 2 in . (13)
With equations (7 -13) we propose a description of the fundamental dynamics of the system. Since we have experimentally estimated α and equation ( 7) relating the initial growth velocity to injection pressure, we can estimate the evolution of the channel length during the invasion for a given injection pressure and final channel length. in and tc = ((α -1)/α) • (x f /v0(Pin)), where C = 10 -3 2 cm/(s•kPa

2 ) and α = -2.5.

IV. DISCUSSION AND CONCLUSION

The channel formations seen in our experiments result in patterns very similar to other natural patterns arising from Laplacian growth, e.g. DLA clusters, viscous fingers in empty Hele-Shaw cells, manganese dendrites, or lightning bolts. In such systems, ramified structures expand at a rate proportional to the gradient of a Laplacian (∇ 2 φ = 0) potential field. This type of patterns is in the DLA universality class, where a fractal dimension of D = 1.71 is expected [START_REF] Løvoll | Growth activity during fingering in a porous hele-shaw cell[END_REF][START_REF] Toussaint | Influence of pore-scale disorder on viscous fingering during drainage[END_REF]. However, the fractal dimensions found for our patterns mainly take values of D ∈ [1.53 -1.60], which is more similar to viscous fingers in porous Hele-Shaw cells [START_REF] Løvoll | Growth activity during fingering in a porous hele-shaw cell[END_REF][START_REF] Toussaint | Influence of pore-scale disorder on viscous fingering during drainage[END_REF]. It has been established that flow in porous media is better described by another Laplacian model, i.e. the Dielectric Breakdown Model (DBM), where the interfacial growth rate is proportional to the pressure gradient of a power η higher than 1, i.e. v ∝ (∇P ) η , where η = 2 for viscous fingers in porous media [START_REF] Løvoll | Growth activity during fingering in a porous hele-shaw cell[END_REF][START_REF] Toussaint | Influence of pore-scale disorder on viscous fingering during drainage[END_REF]. In our system the potential field is the overpressure in the medium, which diffuses into the medium initially (not Laplacian). However, numerical simulations show that the diffusing pressure field quickly approaches the Laplace solution outside finger tips (after 0.4-0.5 s), and in the entire cell on the order of a second [START_REF] Eriksen | Pressure evolution and deformation of confined granular media during pneumatic fracturing[END_REF]. In addition, we observe a feature of the channel growth which is typical for Laplacian growth systems, i.e. there is an active growth zone outside a frozen structure, due to screening of the potential gradient by the most advanced parts of the structure.

The resulting channels in our experiments are differ-ent than channels formed in similar systems having open outer boundaries, such as in [START_REF] Johnsen | Pattern formation during air injection into granular materials confined in a circular hele-shaw cell[END_REF][START_REF] Johnsen | Coupled air/granular flow in a linear hele-shaw cell[END_REF][START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas[END_REF]. In the open systems, the channels are generally much smoother and does not have more than 1-2 branches, while in the closed system we see ramified structures. This is thought to be a consequence of the outer boundary conditions which directly influence how the medium can be deformed. In both systems, the beads are compacted and pushed in front of the growing channel, but in the open system there is also a decompaction front moving inwards from the open outlet, easing the further displacement of beads. On the other hand, in the closed system, after the compaction front hits the outer boundary (which happens between t = 150 -300 ms [START_REF] Eriksen | Pressure evolution and deformation of confined granular media during pneumatic fracturing[END_REF]) further deformation is more difficult due to compaction. Thus, the local configuration of the pressure field on the tips could play a larger role on where the channels grow in this case. In addition, the overpressures imposed for channeling in our experiments (P in = 40-250 kPa) are one order of magnitude larger than the ones for these open systems (P in = 2-20 kPa). Although the cells in [START_REF] Johnsen | Pattern formation during air injection into granular materials confined in a circular hele-shaw cell[END_REF][START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas[END_REF] are smaller than ours, the one in [START_REF] Johnsen | Coupled air/granular flow in a linear hele-shaw cell[END_REF] is of comparable dimensions. The average fluid drag per unit area on the granular medium of height h is given by ∇P • h, and the average threshold for bead displacement is given by [START_REF] Johnsen | Coupled air/granular flow in a linear hele-shaw cell[END_REF]

σ yield = µ(1 -φ)ρgh = µρ s ρgh, ( 14
)
where µ is the Coulomb friction coefficient, ρ s is the solid fraction of the granular medium, ρ is the density of the bead material, and g is the acceleration of gravity. Equa-tion ( 14) is found by considering the weight and friction of the granular medium on the bottom plate while neglecting forces transmitted through solid contacts, i.e. in-plane stress and friction at the top plate, which is reasonable for the initial medium. Then, we have the requirement for bead displacement as

∇P • h ≥ σ yield ⇒ ∇P ≥ µρ s ρg, (15)
which indicates that the pressure gradient necessary to deform the initial granular medium increases with the friction coefficient, solid fraction and the bulk density of beads. In a system where the injected fluid can be considered compressible (like air in our experiment), the overpressure diffuses into the medium with a diffusion constant [START_REF] Eriksen | Pressure evolution and deformation of confined granular media during pneumatic fracturing[END_REF][START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF]

D = κ (1 -ρ s )cµ f = d 2 (1 -ρ s ) 2 180ρ s cµ f , ( 16
)
where κ is the permeability estimated with the Carman-Kozeny expression, d is the bead diameter, c is the fluid compressibility and µ f is the fluid viscosity. By assuming a solid fraction close to 0.5 for the initial medium, we see that the diffusivity of the overpressure is related to the bead size as

D ∼ d 2 cµ f . (17)
Thus, the bead size and cell length influence whether the fluid injection acts as a body force on the medium or more like a surface force on the air-solid interface, i.e. channels do not form with large enough beads or short enough cell because the effective range of the pressure gradient (skin depth) exceeds the system length, such that beads are primarily displaced by the fluid [START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF]. Smaller beads or longer cells lead to a more localized skin depth such that beads are primarily displaced through solid contacts (the zone of mobilized beads eventually exceeds the skin depth) and channels form.

Further, the cell length limits the maximum length of eventual channels formed [START_REF] Eriksen | Pressure evolution and deformation of confined granular media during pneumatic fracturing[END_REF] (longer cells accommodate longer channels). A wider cell leads to an increase in the average number of fingers [START_REF] Niebling | Dynamic aerofracture of dense granular packings[END_REF], but not necessarily an increase in the finger width [START_REF] Johnsen | Coupled air/granular flow in a linear hele-shaw cell[END_REF][START_REF] Johnsen | Decompaction and fluidization of a saturated and confined granular medium by injection of a viscous liquid or gas[END_REF]. However, the channels could become thicker with increased cell width for flow category 4b since a higher number of fingers could merge together. A decrease in cell gap has a similar effect as an increase in cell width, which leads to a higher number of fingers, and an increase in the cell gap above 2 mm usually introduces 3D effects where the fluid tends to flow mostly in the top layer of the granular medium, and to displace grains preferentially there.

For the growth dynamics, the dependency of v 0 on P in , and the exponent α could be system dependent, i.e. the , where Pcat,1 = 9 kPa is the average injection pressure for category 1, where x f = 0. We see that there is a rather scattered (due to randomness in the initial bead configuration from experiment to experiment), but increasing trend in x f for increasing Pin, but also that the final channel length for the highest injection pressures is limited by the closed outer boundary. The fit (Pin -Pcat,1) 0.68 was found from the averaged x f per pressure for Pin ∈ [START_REF] Mazzini | Strike-slip faulting as a trigger mechanism for overpressure release by piercement structures. implications for the lusi mud volcano, indonesia[END_REF]150] parameters could depend on e.g. bead size and confinement. Thus, it is of interest to investigate further how the dynamics are changing in various confined granular media, and see how the Eqs. (7 -13) apply in those cases.

The dependency of x f on P in is also of interest to investigate further, such that if system dependent parameters are known, the expected result of an air injection could be estimated by the injection pressure alone. A brief discussion of what we could find from our experiments is presented in figure 15. Another interesting feature found to be common for all experiments in the system is the typical profile of the finger tips, corresponding to a growth in finger thickness w as function of distance d behind the tip as w ∝ d β where β = 0.68 on average. As shown in figure 11, the fractal dimension of the channels is found to be fairly stable along the main part of the structures, even for the thick fingers in category 4b. The local box-counting method indicates that when there is little or no erosion, as is the case for categories 3a,b and 4a, the final structures end up with typical fractal dimensions between D = 1.53 and 1.60. This range is similar to the one observed for viscous fingers in saturated porous media, i.e. D = 1.53 -1.62, sug-gesting that the channels formed in our experiments are in the DBM universality class of patterns where the interface expands at a rate proportional to (∇P ) 2 . The thicker category 4b structures have a fractal dimension of D = 1.76, which is closer to DLA-patterns and viscous fingering in empty saturated Hele-Shaw cells (D = 1.71). However, the higher fractal dimension is probably also an effect of box-counting inside the main channel which is rather space filling compared to the system size, and this makes the box-counting slopes for category 4b more uncertain than for thinner fingers. Furthermore, the global box-counting dimensions D B are slightly lower than the locally estimated fractal dimensions, which could be a finite size effect since less developed outer parts are also taken into account. However, for the most developed and best preserved channels (Cat. 3b and 4a), we find that D B = 1.53 and 1.54, which is still in good agreement with the fractal dimensions for viscous fingers in porous media. Finally, the fractal dimensions of the front is found to lie between D F = 1.41 and 1.52 for all categories. This indicates that even if the structures themselves appear to change with increasing injection pressure, the roughness of the air-solid interface seems to be more or less the same.

FIG. 1 .

 1 FIG.1. Sketch of the experimental setup. (a) Top-down view of the prepared cell: The granular medium is confined inside the cell by three impermeable boundaries and an airpermeable boundary on the outlet side. The granular medium is placed against the outlet boundary such that it has a linear interface against a region empty of beads on the sealed inlet side, where pressurized air can be injected. (b) Side view of the setup showing the high speed camera placed above. The glass plates are clamped together with aluminum framing while the cell gap is controlled with spacers, which are not shown here.

FIG. 2 .

 2 FIG. 2. Example of image processing. (a): The initial granular medium, (b): A snapshot during the experiment, (c): The resulting binary image where the pattern has been segmented. The segmentation is done by setting all pixels having 30 % less grayscale intensity in a snapshot compared to the initial image to the value 1 while the other pixels are set to the value 0.

  FIG. 3. Typical final structures per category, where the granular medium is gray and areas empty of beads are black, showing that for increasing Pin the fingers become longer and thicker. A full-size figure of the patterns shown here is found in the supplementary material [59].

FIG. 4 .

 4 FIG.[START_REF] Moayed | In-situ stress measurements by hydraulic fracturing method at gotvand dam site, iran[END_REF]. Distribution of new growth over the distance behind the most advanced finger tip, at different intervals, during typical experiments in categories 3a,b and 4a,b. The new growth is measured on intervals between successive snapshots where the channel length grows 10 % of the final channel length for each snapshot. The large initial new growth (0 % -10 %) in all plots corresponds to the initial bulge that forms before the channeling instability begins. The plots indicate the depth of the active growth zone, which goes from very narrow for category 3a (the growth is focused on the most advanced tips), to the whole channel length for category 4b.

FIG. 5 .

 5 FIG.5. Average finger thickness w as function of depth x into the medium for the final structures, averaged per category where channels occur. We see that for the structures in Cat. 4a (erosion) and Cat. 4b (erosion and finger merging) the fingers are thicker -all the way up to the finger tip for Cat. 4b.

FIG. 6 .

 6 FIG. 6. Snapshots from an experiment with Pin = 200 kPa, showing examples of erosion and finger merging. The green rectangle indicates the area cut out and enlarged in the bottom right inset: Here we see erosion inside the channel as beads (gray area) are removed over time. The white arrows point at fingers that will merge into the main channel in the following snapshot, carving out trapped clusters of beads which are then mobilized in the flow direction. A full-size version of this figure is found in the supplementary material [64].

FIG. 7 .

 7 FIG. 7. The average number of fingers N f as function of depthx into the granular medium, averaged over the final patterns in the same category. The plot indicates that for increasing injection pressure, an initial region with fewer fingers extends deeper into the medium before an increase in N f towards the final depth. This suggests more branching of the channels as the beads become harder to displace. The sudden drop in N f close to the tips is due to the finite size of the patterns.

  (b), which shows that for all experiments, β fluctuates around the estimated value β = 0.68 over time, indicating that the tip shapes does not seem to change significantly during the different stages of the experiments.

FIG. 8 .

 8 FIG. 8. (a)The finger thickness w plotted in a log-log plot as function of the distance d from the finger tip. The data shown is the average of about 50 measurements in time for each experiment, averaged again for each category. The data is collapsed by the finger thickness at the base of the tips w b found as the average finger thickness where w has reached a more or less constant level, which seems to be a typical length scale for the tip region. The average slope at different times during several experiments is β = 0.68, indicating a more rounded shape of the tips. In the inset: Snapshot of typical finger tips together with profiles corresponding to β = 0.68 (yellow lines). (b) The slope β averaged per category as function of time during experiments. The plot shows that the slope for the tip profile fluctuates around β = 0.68 over time during the experiments, suggesting that the tips are more rounded than cusp-like at all times.

2 FIG. 9 .

 29 FIG. 9. Average box-counting slopes per category, the lower and upper limits are indicated with the vertical lines and mark the typical thickness of the thinnest fingers (1 cm) up to the cell width (32 cm). In the inset: Box-counting fractal dimensions DB averaged per category, the error bars show 2 standard deviations. Patterns in categories 3b and 4a with DB = 1.53 and DB = 1.54 are within the range of earlier measured fractal dimensions for viscous fingers in saturated porous media (between the dashed lines). The less developed patterns in category 3a is below this range, with DB = 1.41, and the thicker category 4b patterns are close above with DB = 1.63. The solid line indicates DB = 1.71 for Diffusion Limited Aggregation patterns.

2 FIG. 10 .FIG. 11 .

 21011 FIG.10. Average box-counting slopes of the air-solid interface of the final structures per category, the lower and upper limits are indicated with the vertical lines and mark the typical thickness of the thinnest fingers (1 cm) up to the cell width (32 cm). In the inset: Box-counting fractal dimensions DF of the solid-air interface per category. The error bars show 2 standard deviations.

3 2 in

 32 . In figure12(b), the finger tip positions per category is divided by P 3 2

2 FIG. 12 2 in

 2122 FIG.12. (a) The finger tip position x(t) as function of time, averaged for all experiments per category. The curves typically have an initial linear growth with time which crosses over to a decay with time. The patterns grow faster and longer for increasing injection pressure, and the initial linear growth is found to scale proportionally with P 3/2 in as indicated by the dashed lines. (b) The finger tip positions x(t) divided by P 3 2 in , showing that the initial growth fall along the same line. The injection pressures are Pin = 64, 110, 134 and 213 kPa, for category 3a,b and 4a,b respectively.

FIG. 13 .

 13 FIG. 13. (a) Log-log plot of the finger tip velocity v(t) as function of time, averaged per category. The data is collapsed along the ordinate axis by dividing v(t) by P3 2

FIG. 15 .

 15 FIG.[START_REF] Talukder | Pliocene to recent mud diapirism and related mudvolcanoes in the alboran sea(western mediterrenean)[END_REF]. The final channel length x f is plotted as function of Pin -Pcat,1, where Pcat,1 = 9 kPa is the average injection pressure for category 1, where x f = 0. We see that there is a rather scattered (due to randomness in the initial bead configuration from experiment to experiment), but increasing trend in x f for increasing Pin, but also that the final channel length for the highest injection pressures is limited by the closed outer boundary. The fit (Pin -Pcat,1) 0.68 was found from the averaged x f per pressure for Pin ∈ [20, 150] kPa and suggest a power-law relationship. Alternatively, the linear fit 0.46 • (Pin -Pcat,1) fits equally well. The average limiting line xmax,avg was found by averaging x f over the experiments with Pin ≥ 200 kPa, to indicate the influence of the cell boundary. Thus, the results from our experiments indicate that x f increase more or less proportionally with Pin until a limiting length xmax,avg ≈ 59 cm is reached due to the presence of the outer boundary.

  FIG.[START_REF] Talukder | Pliocene to recent mud diapirism and related mudvolcanoes in the alboran sea(western mediterrenean)[END_REF]. The final channel length x f is plotted as function of Pin -Pcat,1, where Pcat,1 = 9 kPa is the average injection pressure for category 1, where x f = 0. We see that there is a rather scattered (due to randomness in the initial bead configuration from experiment to experiment), but increasing trend in x f for increasing Pin, but also that the final channel length for the highest injection pressures is limited by the closed outer boundary. The fit (Pin -Pcat,1) 0.68 was found from the averaged x f per pressure for Pin ∈ [20, 150] kPa and suggest a power-law relationship. Alternatively, the linear fit 0.46 • (Pin -Pcat,1) fits equally well. The average limiting line xmax,avg was found by averaging x f over the experiments with Pin ≥ 200 kPa, to indicate the influence of the cell boundary. Thus, the results from our experiments indicate that x f increase more or less proportionally with Pin until a limiting length xmax,avg ≈ 59 cm is reached due to the presence of the outer boundary.

TABLE I .

 I The flow categories with related pressure values; mean overpressure within two standard deviations Pin ± 2σ, as well as the number of experiments N analyzed per category.

	Category	Pin ± 2σ [kPa]	N
	1	9 ± 5	3
	2	26 ± 18	5
	3a	64 ± 50	5
	3b	110 ± 44	6
	4a	134 ± 43	10
	4b	213 ± 50	4

  FIG. 14. Left: The collapsed finger tip positions x = x(t )/v0 • tc as function of t = t/tc, for individual experiments with various injection pressure. The dotted vertical line indicates the cross-over at t = 1, and the dashed curve shows the proposed function fx(t ). Right: Log-log plot of the collapsed finger tip velocities v = v(t )/v0 as function of t = t/tc for the same individual experiments. The dotted vertical line indicates the cross-over at t = 1, and the dashed curve shows the proposed function fv(t ). Recall that the parameters used for the collapse are given by the injection pressure, as v0 = C • P

	50 kPa	50 kPa
	100 kPa	100 kPa
	150 kPa	150 kPa
	200 kPa	200 kPa
	250 kPa	250 kPa
	f x	(t/t c	)	f v	(t/t c	)
						3
						2

ACKNOWLEDGMENTS

This project has received funding from the European Unions Seventh Framework Programme for research, technological development and demonstration under grant agreement no. 316889, ITN FlowTrans, from the University of Strasbourg through an "Espoirs" grant, from the CNRS-INSU through an "ALEAS" grant, and from the LIA France Norway D-FRACT. This work was partly supported by the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644. We thank Alain Steyer and Miloud Talib for their technical support with experimental equipment. We would also like to thank Emily Brodsky, Daniel Koehn, Piotr Szymczak, Christine Putnis, Ulrich Kelka, Stephen Centrella, Amir Sagy, Monem Ayaz, Anke Lindner, Harold Auradou, Bjørnar Sandnes and Øistein Johnsen for nice and fruitful discussions.


Marcel Moura

Jørgen Knut

Eirik Grude Måløy

Renaud Flekkøy

Toussaint

Eirik Måløy

Flekkøy Grude

Porelab

Renaud Toussaint

Verification of a Dynamic Scaling for the Pair Correlation Function during the Slow Drainage of a Porous Medium

Keywords: numbers: 47.55.-t, 47.15.gp, 64.60.Ht

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

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

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

N (r, t) = r -1 f r D t , (1)
where the dynamic exponent D corresponds to the fractal dimension of the invaded front [START_REF] Mandelbrot | The fractal geometry of nature[END_REF][START_REF] Feder | Fractals[END_REF]. The IP simulations employed in that study confirmed the validity of Eq. (1) numerically and found out additionally that the function f (u) presents the unusual scaling

f (u) ∝ u 1.4 if u 1 u -0.6 if u 1 . (2)
In the current work, we have demonstrated these results experimentally, almost 30 years after the original findings of Furuberg et al. [1]. Fig. 1 shows a diagram of the experimental setup. A quasi-2D porous network is formed by a monolayer of glass beads, with diameters a in the range 1.0mm < a < 1.2mm, randomly placed in the gap of a modified Hele-Shaw cell. The beads are kept in place by a pressurized cushion on the bottom plate of the cell. The porous medium is initially saturated with a viscous liquid (wetting phase), composed of glycerol (80% in weight) and water (20% in weight), having kinematic viscosity ν = 4.25•10 -5 m 2 /s, density ρ = 1.205 g/cm 3 and surface tension γ = 0.064 N/m (see Ref. [START_REF] Moura | [END_REF] for additional details). The outlet channel is connected to a syringe pump, from which liquid can be slowly withdrawn at a constant rate of q = 0.0050 ml/min ≈ 0.14 pore/s, thus assuring that the dynamics happens in the capillary regime, in which capillary forces dominate over viscous ones [START_REF] Lenormand | [END_REF]. Air (non-wetting phase) enters from a width-spanning inlet channel, which is kept open to the atmosphere. Once the capillary pressure (difference in pressure between the non-wetting and wetting phases) is large enough to overcome the threshold associated with a given pore-throat, the invasion of one or more pores happens, and viscous pressure drops are triggered within the medium, thus dissipating energy [3,45]. Pictures are taken every 34 s which leads to an average number of K ≈ 4.7 invaded pores per image (average pore invasion time of t p = 7.2 s). The capillary number is, Ca = ρνq/Σγ = 6.1 • 10 -7 , where Σ = 1.1 • 10 -4 m 2 is the cross section of the model. The choice of a slow constant withdrawal rate reflects the fact that in the IP algorithm employed in Ref. [1], viscous pressure drops are neglected and exactly 1 pore (numerical site) is invaded per time step.

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

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

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

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

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

P t (Θ) ∝ 1 t Θ t -τ all b H (Θ -t) , (3)
where H(x) = 1 for x > 0 and H(x) = 0 otherwise. The dependency on t in the prefactor follows from the normalization in the interval t < Θ < ∞ [47]. We remind that in the IP model, the time corresponds to the size (mass) of the invaded cluster, measured in number of pores. For a cluster of size Θ, the distribution Q Θ (r) of distances r between pores in that cluster is also assumed to follow a power law,

Q Θ (r) ∝ Θ -(1+α)/D r α H Θ 1/D -r , (4)
for r in the interval 0 < r < Θ 1/D , where again the dependency on Θ in the prefactor is obtained by the normalization in the interval 0 < r < Θ 1/D [47]. The pair correlation function N (r, t) is then given by

N (r, t) = ∞ 0 P t (Θ)Q Θ (r)dΘ , (5)
which, using Eqs. ( 3) and ( 4) leads to

N (r, t)r = f r D t =      r D t 1+α D if r D t 1 r D t 1-τ all b if r D t 1 . (6)
The insightful argumentation provided by Roux and Guyon [47] naturally generates the correct functional form of Eqs. ( 1) and ( 2). Nevertheless their analysis for the numerical values of the exponents presents some inconsistent assumptions, as noted by Maslov [48]. The exponent τ all b characterizing the distribution of all backwards avalanches was shown [48] to be given by τ all b = 3 -τ , where τ = 1 + (D e -1/ν) /D is the cluster size distribution exponent (first derived in [49], and verified experimentally in [50]), D and D e are respectively the fractal dimensions of the growing cluster and its boundary and ν = 4/3 is the exponent characterizing the divergence of the correlation length [START_REF] Stauffer | Introduction to percolation theory[END_REF][START_REF] Feder | Fractals[END_REF]. Using the values D = 1.82 and D e = 4/3 (external perimeter growth) [START_REF] Stauffer | Introduction to percolation theory[END_REF][START_REF] Feder | Fractals[END_REF], we find that the short-time exponent is

1-τ all b = -1+ D e -1/ν D =⇒ 1-τ all b = -0.68, (7)
which is consistent with our measurements and very close to the value -0.6 reported in Ref. [1], see Fig. 4. The computation of the long-time exponent brings additional challenges due to the difficulty in estimating the exponent α in Eq. ( 4) [47]. With that in mind, we take an alternative approach here. Consider the situation in which a pore-throat that gives access to a pore at position r 1 has been reached by the liquid air interface at a time t * . From that moment on, that pore-throat and the pore-body at r 1 are available to the invasion. Let us focus here on the case in which this particular pore-throat has a relatively high value of capillary pressure threshold, such that its invasion will typically take a long time to occur. This type of event lies in the long-time regime, for which t r D , and we are interested in addressing their relative probability of occurrence.

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

N (r, t) ∝ 1 (t 1 -t * ) 1 t De/D ≈ 1 t ξ , (8)
where ξ = 1 + D e /D and, in the limit of large times, t 1 -t * ≈ t 1 -t 0 = t. The r dependence in the previous equation is obtained by considering again Eq. ( 1). Since the product N (r, t)r is a function only of the reduced variable r D /t, we have that in the long-time regime,

N (r, t) ∝ 1 r r D t ξ . (9)
Using the literature values D = 1.82 and D e = 4/3 [START_REF] Stauffer | Introduction to percolation theory[END_REF][START_REF] Feder | Fractals[END_REF], we find ξ = 1.73. This value is higher than the exponent ξ = 1.4 reported by Furuberg et al. [1], but is consistent with the experimental data, particularly in the extreme left region of Fig. 4, where t r D and the approximations of the model should hold best.

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

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

FIG. 1 .

 1 FIG. 1. (color online) Diagram of the experimental setup. On top we show an image of the flow and in the detail a zoom in section of the porous network, the defending wetting phase (liquid, blue) and the invading non-wetting phase (air, white).

FIG. 2 .

 2 FIG. 2. (color online) Spatiotemporal map of the invasion up to breakthrough (the average flow direction is from left to right). The colormap shows the elapsed time for the invasion of a given pore (in seconds). The experiment lasts ≈ 33h.

FIG. 3 .

 3 FIG. 3. (color online) Pair correlation function N (r, t) as a function of the distance r for 9 times t (shown in the legend).

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