Germán Varas

Valérie Vidal

Jean-Christophe Géminard

Venting dynamics of an immersed granular layer

Keywords: number(s): 83.80.Fg, 47.57.Gc, 47.85.Dh

Je voudrais d'abord remercier tous les membres du jury pour avoir accepté de participer à cette thèse : à Eric Clément pour avoir présidé la commission, à Yoël Forterre et Philippe Gondret pour le temps dédié au rapport qu'ils ont fait, et à Daniel Schweich qui a été comme un troisième rapporteur. Merci à vous quatre pour votre écoute, vos corrections et vos précieuses opinions qui ont tant apporté lors de la présentation de cette thèse.

Valérie... qui l'aurait dit !! il y a quelques années au Chili j'assistais à une présentation que tu faisais sans savoir ce que l'avenir avait prévu pour moi. Merci pour ton enthousiasme et dévouement, toujours prête à répondre à mes questions avec un sourire. Merci pour ton énergie et l'affection non seulement envers moi mais aussi pour tous ceux qui me sont chers. Pendant ces quatre années tu étais beaucoup plus qu'une directrice de thèse, tu étais une amie.

Jean-Christophe... le 'MacGyver' de la physique ! toujours avec de nouvelles idées pour aborder tous les problèmes d'une manière plus simple, élégante et astucieuse. Tu as toute ma gratitude, admiration et affection.

Un grand remerciement aussi à tous les membres du labo, les secrétaires et les personnes de l'atelier mécanique, spécialement Franck Vittoz pour son aide précieuse avec les différents dispositifs expérimentaux. Aux occupants de mon premier bureau : Nicolas, Sébastien, Vincent et Baptiste pour m'avoir accueilli si chaleureusement tout au début de mon séjour (et aussi pour les après-midis de CS). Un grand merci à l'équipe I, spécialement à Eric Freyssingeas pour son amitié et sa bonne humeur. À ceux qui jour après jour partageaient non seulement la science, mais aussi les petits details du quotidien .... à Osvanny (Papi Ricky !), Ramón, Baptiste et Stéphane... j'espère que ce sera seulement le début ! Enfin, je voudrais remercier mes deux familles 'Varas-Arias' pour tout le soutien durant toutes ces années loin de la maison, et à Fran... qui a été le moteur principal de cette expérience. Maintenant on commence une nouvelle aventure... en famille.

Crater formation

; Lohse et al.

• In Chapter V, we study the dynamics of a crater formation at the free surface of the granular bed. We quantify the size dependence as a function of the flow-rate Φ and the grains size d. Finally we propose a model which accounts for the crater growth.

• In addition to the experiments in immersed granular material, part of this thesis was devoted to another physical phenomenon: the dynamics of bubble bursting.

Parameters

General introduction

Displacement processes can be found in many natural systems and industrial applications. For example, displacement processes in natural porous media or industrial synthetic porous matrices arise in many fields of practical interest, such as petroleum engineering, ground water hydrology, agricultural engineering, soil mechanics and chemical engineering. Among all these examples, it is interesting to analyze in more details two of them that focused our attention from an economical and ecological point of view:

1. Oil extraction: The displacement of one fluid by another is the primary mechanism by which hydrocarbons migrate and are recovered. One of the major problem of oil extraction is that around 40% or more of the porous underground reservoirs is left in a typical oil reservoir by traditional recovery techniques. In order to maximize the extraction from the source rock, gas is directly injected into the reservoir, thus displacing the oil inside. This technique is called "Enhanced Oil Recovery" (EOR) and helps to obtain between 50-60% or more of the reservoir's original oil [START_REF] Lake | Enhanced oil recovery[END_REF].

2. Climate change: Gases are emitted at the sea floor at hydrothermal vents. They are commonly found near volcanically active places and can released huge quantities of carbon dioxide. When the gases are expelled quickly enough, they are thought to potentially induce global climate changes [START_REF] Svensen | Release of methane from a volcanic basin as a mechanism for initial eocene global warming[END_REF]. This mechanism has been proposed as one possible explanation for the climate change in the Toarcian1 . This stage was characterized by a major perturbation of the global carbon cycle expelled by an hydrothermal vent complexe [START_REF] Svensen | Hydrothermal venting of greenhouse gases triggering Early Jurassic global warming[END_REF]. Consequently a better understanding on this process may lead to a better understanding of the causes of the abrupt climate changes [e.g. [START_REF] Wignall | Large igneous provinces and mass extinctions[END_REF]; [START_REF] Courtillot | On the ages of flood basalt events[END_REF]]. Current estimates of global methane fluxes from the seabed to the atmosphere vary between 0.4 and 48 Tg/yr (teragram/year) [START_REF] Judd | Natural seabed gas seeps as sources of atmospheric methane[END_REF].

In addition to the above examples, there are many geophysical phenomena associated with the passage of a fluid through the sea floor [e.g. [START_REF] Mazzini | Strike-slip faulting as a trigger mechanism for overpressure release through piercement structures. Implications for the Lusi mud volcano, Indonesia[END_REF]; Brown, Field, et al. (2008); Porritt et al. (2008aPorritt et al. ( , 2008b))]. One particular application in fluid transport is air sparging. In this process, air is injected into the subsurface below the lowest known depth of contamination and due to buoyancy (and other mechanisms such as volatilization and biodegradation), the air helps to remove or degrade the contaminants [e.g. [START_REF] Semer | An experimental investigation of air flow patterns in saturated soils during air sparging[END_REF]; [START_REF] Peterson | Airflow geometry in air sparging of fine-grained sands[END_REF]; [START_REF] Reddy | Effects of soil heterogeneity on airflow patterns and hydrocarbon removal during in situ air sparging[END_REF]; [START_REF] Adams | Laboratory study of air sparging of TCE-Contaminated saturated soils and ground water[END_REF]].

The observations presented here show a complex behavior which is not easy to describe. For example, acoustic sea-floor backscattering has shown that when there is an homogeneous reservoir depth below the bed sea, the gas (or fluid) emission at the surface can be concentrated in singular points or more interestingly exhibit clusters formation (Naudts et al., 2008) (Fig. 1.1). It would be interesting to invert the problem and infer from the observation at the surface, some information about the characteristic of the reservoir at depth.

In this introductory chapter we present some examples of geological structures related to the displacement processes described above (Sec. 1.1.1). As all the phenomena presented involve a complex medium, in particular a pore structure, we then describe in details previous studies of fluid flow in porous medium (Sec. 1.2). In particular we describe previous experiments based on the injection of a fluid in a granular material immersed in the same fluid (two phase flow), when the injection point is punctual at the center or homogenous along the base (Sec. 1.3.1). We then present the phase diagram for a three phase flow configuration, and fully describe the dependence of the invasion pattern on the parameters (Sec. 1.3.2). Finally, we describe the outline of this PhD thesis.

Piercement structures

When a pressurized fluid is forced through a rock or sediments it generates geological patterns that are called "piercement structures". They exist in a range of different forms and are observed worldwide. Here, we describe briefly the formation and mechanism of four different piercement structures that can be observed in nature [Fig. 1.2]:

• Pockmarks : Pockmarks are crater-like depressions in the seabed [Fig. 1.2 (a)].

They are broad, circular depressions typically 30 to 40 m across by 2 or 3 m deep.

The formation is generally assumed to be caused by the expulsion of fluids from the sediment layer [START_REF] Kvenvolden | Seabed pockmarks and seepages -impact on geology, biology and the marine-environment[END_REF]. The agent is supposed to be of various gases, especially methane, originating from depth [START_REF] Kelley | Giant sea-bed Pockmarks -Evidence for gas escape from Belfast bay, Maine[END_REF], but may also be fresh groundwater [START_REF] Khandriche | Freshwater induced pockmarks in Bay of Eckernfoerde, Western Baltic[END_REF] or pore water expulsed by compaction [START_REF] Harrington | Formation of pockmarks by pore-water escape[END_REF]. They have been found in a number of widely spread localities [e.g. [START_REF] Hovland | Complex pockmarks with carbonate-ridges off mid-Norway: Products of sediment degassing[END_REF]; [START_REF] Newman | Active methane venting observed at giant pockmarks along the us mid-atlantic shelf break[END_REF]; [START_REF] Hammer | Piston coring of Inner Oslofjord pockmarks, Norway: constraints on age and mechanism[END_REF]].

• Hydrothermal vent complexes: The gas leakage from a surface vent or the fluidization of initially consolidated sediments are representation of hydrothermal vent complexes [Fig. 1.2 (b)]. Some examples are found in the Karoo Basin, South Africa [START_REF] Svensen | Structure and evolution of hydrothermal vent complexes in the karoo basin, south africa[END_REF]; [START_REF] Lock | Discussion on structure and evolution of hydrothermal vent complexes in the Karoo Basin, South Africa[END_REF]] and the Norwegian sea [START_REF] Svensen | Seep carbonate formation controlled by hydrothermal vent complexes: a case study from the Voring Basin, the Norwegian Sea[END_REF]. As mentioned before, the comprehension of this phenomena could help to understand the carbon cycle perturbation [START_REF] Svensen | Release of methane from a volcanic basin as a mechanism for initial eocene global warming[END_REF]).

• Mud volcanoes : They are the geologically manifestations of vertical fluid flow and mud eruption in sedimentary basins [Fig. 1.2 (c)]. Their formation is the result of the expulsion of an overpressurised complex mixture of sediments, water, and various chemicals where hydrocarbon gases (mainly methane) and in some cases petroleum are the dominant components. Thus the presence of seepage indicates the existence of subsurface accumulations of hydrocarbons and migration pathways. They are present worldwide in different zones as in the Egypt offshore [START_REF] Mastalerz | Active venting at the Isis mud volcano, offshore Egypt: Origin and migration of hydrocarbons[END_REF], Azerbaijan [START_REF] Planke | Mud and fluid migration in active mud volcanoes in Azerbaijan[END_REF] and Indonesia [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 through piercement structures. Implications for the Lusi mud volcano, Indonesia[END_REF].

• Kimberlite pipes : Kimberlite pipes are generally formed from fragmentation of deep dyke complexes [Fig. 1.2 (d)] [START_REF] Lorenz | Maars and diatremes of phreatomagmatic origin: A review[END_REF][START_REF] Clement | The origin of kimberlite pipes: An interpretation based on synthesis of geological features displayed by southern african occurrences[END_REF]. This last mechanism may also explain the formation of the phreatomagmatic complexes in the Karoo Basin [START_REF] Surtees | The evolution of the Brosterlea Volcanic Complex, Eastern Cape, South Africa[END_REF]; [START_REF] Mcclintock | The volcaniclastic opening phase of Karoo flood basalt volcanism, Drakensberg Formation, South Africa[END_REF]]. Numerous studies have been devoted to this phenomenon (Porritt et al., 2008a;Brown, Field, et al., 2008). Kimberlites pipes are the most important source of primary diamonds [START_REF] Walters | The role of fluidisation in the formation of volcaniclastic kimberlite: Grain size observations and experimental investigation[END_REF]Gernon, Sparks, & Field, 2008;[START_REF] Torsvik | Diamonds sampled by plumes from the core-mantle boundary[END_REF]. Examples of such phenomenon can be found in Botswana (Brown, Gernon, et al., 2008) and the Southern Africa (Gernon, Sparks, & Field, 2008).

Each one of the above cases describe different geological structures that are produced by the same mechanism: the flow of a fluid through a complex deformable Figure 1.2: Schematic drawing of different piercement structures observed in nature (from [START_REF] Nermoen | Experimental and analytic modeling of piercement structures[END_REF]) -(a) Submarine pockmarks off-shore Norway (after [START_REF] Hovland | Petroleum Geochemistry in Exploration of the Norwegian Shelf[END_REF]. (b) Witkop II hydrothermal vent complex, Karoo Basin, South Africa (after [START_REF] Svensen | Structure and evolution of hydrothermal vent complexes in the karoo basin, south africa[END_REF]. (c) Mud volcano, Azerbaijan (after [START_REF] Planke | Mud and fluid migration in active mud volcanoes in Azerbaijan[END_REF]. (d) Kimberlite pipe, Botswana (after [START_REF] Walters | The role of fluidisation in the formation of volcaniclastic kimberlite: Grain size observations and experimental investigation[END_REF]. matrix. Even if the media can be different (gas emission, sediments, rock fragments, etc) and also, act at different length scale (∼ 100 m for pockmarks and the hydrothermal vent complexes, and ∼ 1 km for mud volcanoes and kimberlite pipes) the geological structures display a similar geometry (pipe formation) and, in particular, due to the emission of material, the formation of a fluidized zone in the 'bulk' and a crater at the surface. (2004); [START_REF] Zheng | Impact craters in loose granular media[END_REF]; [START_REF] Wada | Numerical simulation of impact cratering on granular material[END_REF]; [START_REF] Vet | Shape of impact craters in granular media[END_REF]], as well in immersed granular materials [e.g. [START_REF] Gostiaux | Dynamics of a gas bubble rising through a thin immersed layer of granular material: an experimental study[END_REF]; [START_REF] Rigord | Localized instability of a granular layer submitted to an ascending liquid flow[END_REF]; Zoueshtiagh and Merlen (2007)]. In the latter case, a crater can be generated by underwater rings [START_REF] Suzuki | Modes of ejecta emplacement at martian craters from laboratory experiments of an expanding vortex ring interacting with a particle layer[END_REF] involving fluidized ejecta dynamics, or underwater impact craters generated by landslide [START_REF] Fritz | Landslide generated impulse waves. 1. Instantaneous flow fields[END_REF]. Craters in immersed granular materials can result either from two-phase or three-phase flows. In particular, water or gas flowing through an immersed granular bed can induce localized instabilities and fluidization [START_REF] Rigord | Localized instability of a granular layer submitted to an ascending liquid flow[END_REF]Zoueshtiagh & Merlen, 2007), which eventually lead to the formation of craters at the free surface.

Introduction

In all of these phenomena (piercement structure and crater formation) we have to deal with the complex pore structure of the medium and how it affects the distribution, flow, the displacement of one or more fluids, or the dispersion (i.e., mixing) of one fluid in another. Each process is, in itself, very complex. For example, the displacement of one fluid by another can be carried out by many different mechanisms, which may involve heat and mass transfer, thermodynamic phase behavior and phase change, and the interplay of various forces such as viscosity, buoyancy, and capillarity. If the solid matrix of the porous medium is deformable, its porous structure may change during the flow or any other phenomenon, the system becoming even more complex. Therefore, it is extremely important to understand the physics behind at the pore scale. As physicists we are interested in the individual and global mechanisms behind this phenomenon. In the next chapter we present a brief review of previous works in porous media, upon which further analysis and interpretation will be based.

Porous medium 1.2.1 Description

A porous medium consists of a matrix with a large amount of microscopic pores and throats which are typically narrow tubes where fluid can pass through. In nature we can consider a porous medium as the portion of rock not occupied by solid matter, that is, the void spaces (e.g. pores, interstices and fissures). Even if the global description of the medium is simple, the description of a flow in the system is not simple at all, due to its complexity (usually a network). In principle the description of the flow in a single tube is given by simple equations, but the full description of the tube network is most of the time impossible. The porous medium (matrix) can be static [Fig. 1.3 (a)], i.e. the morphology of the medium does not change during a given process, or deformable [Fig. 1.3 (b)], i.e. its morphology changes due to chemical reactions or to physicochemical interactions between the pore surface and a fluid.

We are interested in the description of a fluid that flows through a porous medium. In a consolidated medium (rigid), the Darcy law relates the permeability and the is created by molding techniques using transparent polyester and photographically etched mould, the distance between two nodes is about 0.8 mm (Lenormand & Zarcone, 1985). (b) Soft porous network. Example of a fluid foam, the surface of a soap froth [START_REF] Graner | Two-dimensional fluid foams at equilibrium[END_REF].

pressure gradient:

q = - k µ ∇P , (1.1)
where q is the flux (discharge per unit area in m/s), k the permeability (m 2 ), µ the viscosity (Pa • s) and ∇P is the pressure gradient vector (Pa/m). In addition to Darcy law, the description of the system at pore scale is completed by the capillary (or Laplace) pressure, which is defined by the difference between the pressures in the non-wetting and wetting phase (Fig. 1.4):

∆P c = P nw -P w = γ 1 R 1 - 1 R 2 (1.2)
where γ is the surface tension between the liquids and R 1 and R 2 are the two main radii of curvatures for the interface.

When capillary forces are predominant over viscous forces, the injection of an immiscible fluid into a porous medium filled with another fluid produces two types of displacement process depending on the wettability of the fluids (Fig. 1.5):

Wetting fluid Nonwetting fluid P w P nw Figure 1.4: Wetting phases at pore scale (from [START_REF] Lovoll | Competition of gravity, capillary and viscous forces during drainage in a two-dimensional porous medium, a pore scale study[END_REF] -The difference between the pressure in the nonwetting fluid P nw and that in the wetting fluid P w is given by the Laplace law; during drainage, capillary forces act against the displacement, and larger pores are more easily invaded.

Introduction

Drainage Imbibition Depending on the property of the invading and host fluid [e.g. wetting (represented in blue) or non-wetting (represented in white)] we can have two types of invasion. When the non-wetting fluid displaced the wetting, the process is called drainage. The inverse process is called imbibition.

• Imbibition: Occurs when a wetting fluid displaces a non-wetting fluid2 .

• Drainage: Occurs when a non-wetting fluid displaces a wetting fluid.

The relative wettabilities, viscosities, and densities of the fluids, as well as the heterogeneity of the underlying porous media, play an important role in the competition process. In both cases (imbibition or drainage), the displacement process generates patterns formation between the interface of the two fluids, and the different structures obtained can be divided into three major flow regimes: viscous fingering, stable displacement and capillary fingering (invasion percolation) [START_REF] Sahimi | Flow phenomena in rocks -from continuum models to fractals, percolation, cellular-automata, and simulated annealing[END_REF]. The mechanisms of the displacements in drainage and imbibition are quite different and the two cases should not be confused [START_REF] Lenormand | Mechanisms of the displacement of one fluid by another in a network of capillary ducts[END_REF]. Typically, slow drainage is characterized by a piston-like motion inside the pores where the invading non-wetting fluid only enters a pore if the capillary pressure is equal to or greater than the threshold pressure of that pore. The threshold pressure corresponds to the capillary pressure in the narrowest part of the pore. However, in imbibition at low injection rate the invading fluid will enter the narrowest pore before any other is considered.

Dimensionless numbers

A set of dimensionless numbers is usually defined to quantify the relative magnitudes. The Capillary number Ca is the typical ratio of the viscous pressure drop at pore scale to the capillary pressure, while the Bond number Bo quantifies that of the typical hydrostatic pressure drop over a pore to the capillary pressure:

Ca = ∆p visc ∆p cap = η v γ , (1.3)
and .4) where η is the fluid viscosity, v the characteristic fluid velocity, γ the fluid interface tension, ∆ρ the fluid density difference, g the acceleration of gravity, and a the typical pore size.

Bo = ∆p grav ∆p cap = g∆ρa 2 γ , ( 1
During slow displacement, capillary forces dominate over viscous forces leading to patterns described by the invasion percolation model [START_REF] Wilkinson | Invasion percolation -a new form of percolation theory[END_REF]. However, in most cases there are density differences so that the gravity causes hydrostatic pressure gradients [START_REF] Wilkinson | Percolation model of immiscible displacement in the presence of buoyancy forces[END_REF]. These gradients compete with capillary forces that are randomly distributed due to the randomness of the porous medium.

Experimental and numerical studies in porous media

The physical description of gas flow injection into saturated rigid porous media has focused attention over the last decade. Many experimental works have been performed on the migration of air channels [START_REF] Kong | Migration of air channels: An instability of air flow in mobile saturated porous media[END_REF], the morphodynamics [START_REF] Kong | Morphodynamics during air injection into water-saturated movable spherical granulates[END_REF] and the local geometry of gas injection into saturated homogeneous porous media [START_REF] Selker | The local geometry of gas injection into saturated homogeneous porous media[END_REF]. The transition from incoherent to coherent buoyancydriven system [START_REF] Geistlinger | Direct gas injection into saturated glass beads: Transition from incoherent to coherent gas flow pattern[END_REF], or a system destabilized by gravity have also been investigated [START_REF] Glass | Gravity-destabilized nonwetting phase invasion in macroheterogeneous porous media: Experimental observations of invasion dynamics and scale analysis[END_REF]. The modeling at pore-scale of gas flow pattern [START_REF] Geistlinger | Pore-scale and continuum modeling of gas flow pattern obtained by high-resolution optical bench-scale experiments[END_REF] as well as size-dependent scaling of capillary invasion [START_REF] Hirsch | Size-dependent scaling of capillary invasion including buoyancy and pore-size distribution effects[END_REF] has been studied. Pattern formation is also observed during slow drainage of a granular-fluid in two-dimensional confinement [START_REF] Sandnes | Labyrinth patterns in confined granular-fluid systems[END_REF]. At the pore-network level, this problem can be modeled by invasion percolation (IP) [START_REF] De Gennes | Lois générales pour l'injection d'un fluide dans un milieu poreux aléatoire[END_REF][START_REF] Chandler | Capillary displacement and percolation in porous-media[END_REF], in which the front separating the two fluids advances by penetrating the pore throat at the front with the largest size (smallest capillary resistance). The properties of IP and its close connection to ordinary percolation (OP) have been extensively studied.

As the global process is difficult to explain, numerical simulations were extensively used to mimic the system. The quasi-static gravity-destabilized displacement process can be simulated using algorithmic models based on the invasion percolation model. Among these we can find the Diffusion-limited-aggregation (DLA) (T. Witten & Sander, 1981), Eden model (Eden, 1961), Vold-Sutherland model (Vold, 1963;Sutherland, 1966), etc. These models generate structures similarly to the fluid-fluid invasion process in a 2D porous media (Fig. 1.6). Percolation theory in numerical simulation is used to describe the displacement of one fluid by another when capillary forces are important compared to viscous forces. In imbibition, the capillary domain is described either by a compact cluster growth (small aspect ratio) or percolation theory (large aspect ratio). For example, fast viscous-dominated displacement in the limit of large viscosity ratio (viscous fingering) is modeled by diffusion limited aggregation (T. A. [START_REF] Witten | Diffusion-limited aggregation[END_REF][START_REF] Patterson | Diffusion-limited aggregation and two-fluid displacements in porous media[END_REF]. Invasion percolation (IP) [START_REF] Wilkinson | Invasion percolation -a new form of percolation theory[END_REF] describes the opposite extreme of slow displacement, where capillary forces dominate over viscous and lead to percolation patterns.

Finally, a possible way to differentiate the processes is to quantify the fractal dimension in invasion percolation (Lenormand & Zarcone, 1985). Theoretical argument can be developed for the fractal dimension calculations in OP [START_REF] Feder | Fractals[END_REF], leading for two dimensional system to D = 91/48 ≃ 1.89. When a trapping rule is applied, regions surrounded by the invading phase become disconnected, and experiments and simulations have reported that the fractal dimension in IP is D ≃ 1.82 [START_REF] Wilkinson | Invasion percolation -a new form of percolation theory[END_REF]Lenormand & Zarcone, 1985).

Until now we have described the most basic possible situation, a rigid porous medium. But natural phenomena present unconstrained media which are more complex due to the displacement of one or more fluids and the deformable matrix. In the next section we present several works which aim at describing the dynamics in unconstrained media.

A C B

Unconstrained medium 1.3.1 Two phase flows

A more realistic approach to understand complex natural systems (Sec. 1.1), usually unconstrained and with a free surface, consists in modeling them with cohesionless material which are easier to understand. A typical experiment consists in filling a cell with a granular material immersed in a fluid (note that we consider fluid in a general form, that is, liquid or gas). The same or another fluid is then injected at the bottom of the recipient through multiple or single injection points. In the following, we synthesize the dynamics observed in some examples for the case of single or multiple injection points for a two phase system. Several experiments of fluid injection through multiple points at the bottom of a container filled with the same fluid and sand mixtures have been performed [START_REF] Wilhelm | On the onset of flow instabilities in granular media due to porosity inhomogeneities[END_REF][START_REF] Mörz | An experimental setup for fluid venting in unconsolidated sediments: New insights to fluid mechanics and structures[END_REF][START_REF] Nermoen | Morphological transitions in partially gas-fluidized granular mixtures[END_REF]. In all cases, the results share the same three main flow regimes which depend on the flow-rate and the hydraulic gradient 3 (Fig. 1.7). At low flow-rate the grains form a fixed skeleton through which the fluid percolates homogeneously in an upward direction [Fig. 1.7 (a)]. The porosity distribution is the same within the whole region. Increasing the flow-rate, larger pores and small conduits with lengths much shorter than the height of the sample form at the soil surface [Fig. 1.7 (b)]. For high flow-rate, larger pipes begin to form, crossing the soil from bottom to top [Fig. 1.7 (c)]. Within the pipes, grains and water move like a fluid. At the surface, the pipes erupt in small volcanoes. These pipes are not stationary, and migrate across the surface, progressively fluidizing more and more the surrounding regions. Figure 1.8 shows an example for the particular case where the solid phase consists in a bimodal granular mixture and the fluid is air [START_REF] Nermoen | Morphological transitions in partially gas-fluidized granular mixtures[END_REF]). We observe, in a similar way, the three main flow regimes which depend on the injection flow-rate 4 .

A C B

Finally the same regimes were also reported in the case of a punctual injection of water in a water-saturated granular bed (Zoueshtiagh & Merlen, 2007). The injection of water is located at the bottom of the cylinder. Figure 1.9 shows one example of fluid emission at the granular bed surface. The photograph displays sand pattern developing at the sand-water interface as a result of water flowing vertically from beneath. Note that even if the injection is punctual, the gas emission at the surface is located over a broad region.

Flow regimes

It is interesting to note that even if the experimental setups presented above (3D and 2D) have different injection types and phases, the global behavior of the system is the same. The dynamics of an unconstrained saturated granular medium crossed by the same fluid can thus be summarized in the three following regimes:

1. Darcy flow: At low flow-rate, the fluid crosses the granular bed as if it were a rigid porous medium.

2. Formation of microchannels (or conduits): Increasing the flow-rate, small heterogeneities are created. Small conduits appears and some events are observed at the surface.

3. Hydraulic failure: Above a critical flow-rate, the system is fluidized and a channel is formed. At the surface, small eruptions are observed.

Three phase flows

When adding one more phase to the system, we get closer to the real scenario. In this section we describe two works which report a tentative phase diagram on this type of system.

Invasion pattern phase diagram

Recently, [START_REF] Sandnes | Patterns and flow in frictional fluid dynamics[END_REF] presented the study of a generic granular material that settles in a Newtonian fluid. They slowly inject air into a linear Hele-Shaw cell5 loaded with polydisperse glass beads (∼ 100 µm diameter) immersed in a water/glycerol solution. They explored the displacement dynamics of the system and mapped the dynamic response of the frictional fluid as it is displaced by the air, focusing on the effects of the granular filling fraction ϕ, the system stiffness K and the injection rate Φ. Several new dynamic modes were uncovered. Fig. 1.10 displays the [Φ, ϕ -1 ] phase space diagram. Depending on the parameters, different patterns are observed. At low filling fraction the injected air advances in a slow and quasi-continuous fingering process where side-branching of fingers and a random growth direction produces a branched labyrinthine structure. At low flow-rate viscous forces are negligible, thus capillary and friction forces are thereby effective. For high filling fractions, the system changes from a fluid to a solid-like behaviour. The transition to fracturing is observed from approximately ϕ > 0.9, where the poorly compacted granular material can be described as a deformable, or weak, porous medium. When the system approaches the close packing limit, the granular material resembles a rigid porous medium, and the displacement of the fluid progresses pore by pore, governed by local capillary and viscous forces.

A similar studied was performed by [START_REF] Chevalier | Morphodynamics during air injection into a confined granular suspension[END_REF]. In this case, the morphodynamics during air injection into a confined granular suspension6 was studied when a less viscous fluid (air) pushes a more viscous fluid in a thin linear channel (or

Introduction

Hele-Shaw cell). When increasing the volume fraction ϕ a transition occurs (Fig. 1.11). For lower volume fraction ϕ, we observe a Saffman-Taylor instability. Increasing the grain fraction the resulting pattern formed by the air transits from: (a) the classical Saffman-Taylor instability, to (b-c) pattern formation in suspension and finally (d) viscous fingering in a classical porous medium.

Investigation proposal

In all the examples presented here, the study of the gas emission was generally focused at short time scale. In general the invasion pattern is characterized as a function of the flow-rate and the compaction. In particular we aim at describing the dynamics of air crossing an unconstrained immersed granular medium, and how the gas emission will be located at the free surface7 , when injecting from a punctual inlet at the bottom of the cell. The system is studied experimentally, numerically and theoretically.

Granular material

Granular materials are the perfect choice to model in the laboratory this type of phenomena. They provide a well-controlled and reproducible material to work on. In order to characterize the granular material used in our experiments, we analyzed several samples of glass beads (USF Matrasur, sodosilicated glass) which are previously sieved in order to control their size (diameter d). We determine the grain size distribution for each batch (given by the provider) by means of a microscope (Leitz, Laborlux). We take several pictures from different samples and then calculate with the assistance of a computer program (Matlab, MathWorks) the size distribution for a large number of grains (typically around 500 grains). The results indicate that the three grain size d =200-250, 250-425, 400-500 µm are well-described by a Gaussian distribution 8 (Fig. 1.12). The values for the average µ and the correspondent standard deviation σ for the grains size are summarized in table 1.4.1. We define the average grain size as d = µ ± σ for each series of grains. Note that in the following chapters, we used the batch name to refer to the grain size d.

Outline

In this thesis we present the experimental study of the dynamics of air crossing an unconstrained immersed granular medium in a tank (3D) and a Hele-Shaw cell (2D), as well a theoretical and numerical modeling for this system. The different points which we develop in this thesis are the following: 433.9 ± 1, 5 40.1 ± 3.5

• In Chapter II, we study the gas emission location at the surface in a 3D configuration when injecting air at constant flow-rate. Also, we describe the geometry produced by the successive generation of branches in a 2D configuration. We also develop a simple analytical model and contrast the results with numerical simulations.

• In Chapter III, we extend the numerical simulations, analyzing the morphology of the invasion pattern for a constant volume of air injected. We compare the results for two types of boundary conditions.

• In Chapter IV, we present preliminary experimental results of the evolution of the fluidized zone contour. We show the evolution of the diffusion coefficient as a function of the tilting angle θ and calculate the velocity of the front.

In this chapter we present the experimental study of the gas emission location at the surface of a granular bed immersed under water in a 3D configuration, and the geometry formed by the successive generation of branches in a 2D configuration. After describing the 3D and 2D setups (Sec. 2.1), we analyze the emission of bubbles at the surface (Sec. 2.2). We then analyze the invasion regimes observed at a short time scale in the 2D setup (Sec. 2.3.1) and the geometry formed at long time scale (Sec. 2.3.2). We develop a simple analytical model (Sec. 2.4) and contrast the results with numerical simulations (Sec. 2.5). Finally, we discuss the generation of clusters in the 3D (Sec. 2.6).

Experimental setups

3D experimental setup

The experiment consists in injecting gently air at the base of a column of immersed grains and observe the loci of the resulting bubble emissions at the free surface. The setup is made of a Plexiglas cylindrical tank (24 cm diameter). Air is injected into a chamber connected at the bottom of the cylinder through an inlet hole (Fig 2 .1)1 . The air injection is performed by means of a reducing valve and a capillary tube. The relationship between the air flow-rate and the pressure drop in the capillary tube can be found by the Poiseuille equation:

Φ = πr 4 8ηL ∆P (2.1)
where ∆P is the pressure drop, η the air dynamic viscosity, Φ the flow rate and L (resp. r) the length (resp. radius) of the capillary tube. To ensure that the pressure drop in the capillary tube is much larger than the fluctuations of the pressure inside the chamber (providing that the flow rate is not significantly altered by the dynamics of the granular layer and remains constant) we choose a capillary tube of L = 5 cm and r = 200 µm. These values give us a pressure drop of ∆P ∼ 3 kPa which is sufficiently larger than the typical variation of pressure in the system (Fig. 2.

2).

The flow rate Φ is tuned by changing either the capillary tube (rough) or the pressure difference imposed by the reducing valve (fine) and is subsequently measured in the range of 0 to 5 mL/s (within 0.1 mL/s) by means of a homemade flow meter. A sensor (MKS Instruments, 223 BD-00010 AB), connected to an acquisition board (National Instruments, PCI-6251), monitors the pressure difference, ∆P (to within 10 Pa), between the chamber and the outside atmosphere.

Protocol

The cylinder is filled with the glass beads. The height of the granular bed h g ranges from 2 to 24 cm. Distilled water is added and the material is vigorously stirred with a paddle to eliminate the gas bubbles trapped inside. The free surface of the granular bed is leveled and the water height above is h w ∼ 1 cm, which is small enough to prevent the rising bubbles from hiding the emission locus.

Air is initially injected at a given flow rate Φ. After penetrating the bulk, the air passes through the system and we observe at the surface a series of bubble emissions. The loci of the gas emission (which can vary through time) is detected by a webcam (Logitech, QuickCam S7500) connected to a PC. The cylinder is illuminated with a lateral source (transparency flat viewer, Just NormLicht, Classic Line) that is placed slightly out-of-axis to avoid direct reflexions from the liquid surface into the camera. This type of illumination provides a good contrast for the emission of a bubble (Fig. 2.1b), which can be later analyzed to find its exact position. A free software 0 0.5 1 2.6 (Astra Image Webcam Video Grabber) is used to record videos at 20 images per second with a resolution of 640 × 480 pixels. An automatic detection (Matlab, MathWorks) makes it possible to determine the loci of the gas emission with an accuracy of about 2 mm.

2D experimental setup

The 3D setup does not allow us to visualize the air path inside the granular bed. Therefore, we built another experimental setup which aims at characterizing the dynamics of the air path. The 2D setup consists of a vertical Hele-Shaw cell, with a gap e between the two glass plates (Fig. 2.3a). Similarly to the other setup, we fill the cell with grains (height h g ) immersed in water (height h w above the granular bed). In order to control the gravity effects, the cell can be tilted by an angle θ with respect to the vertical so as to produce an effective gravity g eff ≡ g cos θ.

In order to account for the dynamics of the system, we image from the side. In a first configuration (Chap. 2 and 3), a transparency flat viewer (Just NormLicht, Classic Line) located behind the cell is used to achieve a homogeneous lighting of the whole bed which is imaged with the help of a Webcam (Logitech, QuickCam Express) connected to a PC. A small software (Astra Image Webcam Video Grabber) makes it possible to take one image of the system every 10 seconds and, thus, to record the dynamics during several hours (typically 24 hours).

Air injection

The injection of air through an inlet (2mm diameter) located at the center of the lower edge is insured by a mass-flow controller (Bronkhorst,.

The flow rate Φ can be tuned in the range 0.2 to 2.2 mL/s. For some special cases we used a fixed air volume injection. This configuration is obtained by connecting directly the inlet to a syringe pump (Harvard apparatus syringe Pump 11 Plus) which can be programmed to inject a specific air volume at a given air flow-rate.

Protocol

Initially, the air flow-rate Φ is set to zero. The initial condition in the system consists as follows. First, we pour the grains roughly up to the height h g . We then fill completely the cell with water and close the upper aperture. We then flip the cell bottom-up, which forces most of the small bubbles trapped in the cell to separate from the surrounding grains. While the grains are not fully deposited (very loose packing), we flip the cell back to its original position, therefore ensuring that all the small bubbles are completely extracted from the surrounded grains. We then fix precisely the water (h w ) and grains height (h g ) and level the granular layer with a small paddle. Note that the global packing is controlled but not the local packing.

After opening the valve, we observe that the air is creating one or several paths between the injection nozzle at the bottom and the free surface of the granular bed [Fig. 2.3(b)]. During this transient regime, some air channels merge, or cannot reach the upper layer of grains (in this case, air bubbles remain trapped within the granular layer). After several minutes, the air crosses the granular bed along the vertical and bubbles are emitted in the water.

Gas emission location at the surface (3D)

In this experiment we aim at characterizing the loci of the gas emission at the free surface of the immersed granular bed. We measure the position of the bubble emission for 10 minutes at constant air flow Φ for different heights of grains h g . We thus obtain the position (x, y) at the granular bed surface for the successive bubble emissions (see Fig. 2.4).

General observations

We repeat the experiments for three different grain sizes (d =150-250, 250-425 and 425-600 µm) using three different flow-rates (Φ =1, 2, 3 mL/s). Qualitatively, we observe the following behaviors for the different grain sizes.

Small grain size (150-250 µm)

• When the air is injected, a large air pocket forms below the surface. Typically -at any given air flow Φ -the air pocket grows for 3 to 5 minutes until the air trapped is suddenly liberated. In this case, the bubble emission at the surface occurs over a region limited by the size of the underlying air pocket (a few centimeters). We can understand this behavior as follows: as the dynamics is controlled by the Laplace overpressure γ/d, when the grains are small, the Laplace overpressure is large enough to lift the whole granular column above as it exceeds the typical hydrostatic pressure ρ g g h g , where ρ g stands for the density of the grains. Large grain size (425-600 µm)

• In this range of grain size, the air crosses the immersed grain layer easily, producing not a single, but multiple simultaneous bubble emissions at the surface. This behavior lasts for several seconds, up to several minutes, until stabilizing to a single bubble emission location at the vertical of the injection point. Contrary to the case of the small grain size, the Laplace overpressure is small due to the large grain size and the air flows easily between the grains. The air then invades the material almost in the whole available free space between the grains as a consequence of the system polydispersity.

Intermediate grain size (250-425 µm)

• Bubbles are emitted at the surface of the granular layer, without any largescale surface deformation. The emission consists of successive, single bubbles, apparently randomly distributed in space around the centered air injection point.

In most cases the emissions are not time dependent (contrary to the two other grain sizes) at any flow-rate (Fig. 2.4).

In this work we focus on the intermediate grain size. Indeed, they exhibits venting dynamics similar to the natural phenomena described in chapter 1.

Gas emission distribution

During the experiment, we observe several hundreds of events (bubbles emission) at the free surface (Fig. 2.4). From the image analysis, we obtain the distance r between the bubble emission point and the vertical of the injection point at the surface of the granular layer [Fig. 2.5(a)]. The figure shows the evolution of r through time for a given set of parameters (d, h g , Φ). We observe that the bubble emission location is widely distributed and rather homogeneously distributed around the center. Figure 2.5(c) shows the autocorrelation function R of the distance r. As we can observe, R drops down quickly after a few events, and vanishes after about 10 bubble emissions. This results confirm that there is no long-term correlation, and we can approximate that each bubble emission point is independent from the previous bubble location.

In order to quantify the typical size of the events distribution, we determined the cumulative probability defined as:

F ≡ r 0 p(u)du, (2.2)
where p(r) is the probability to find a bubble from a distance r of the center. We find that the complementary cumulative distribution function 1 -F is rather well described by a Gaussian function [Fig. Reporting ζ for the different grains height h g (Fig. 2.6), we find that the data are compatible with following relationship:

ζ ∝ h g , (2.3)
which does not depend significantly on the air injection flow rate Φ. The characteristic size ζ scales as the square root of the grains height h g for the different air injection flow rate Φ. To investigate the mechanisms leading to this relationship, and observe in details the air path through the granular layer, we performed a series of experiments in 2D which we describe in the next section.

Air path through the granular layer (2D)

The 2D experimental setup allows the directly visualization of the air crossing the granular bed. The observations indicate that the path generated by the air is dynamic and changes in time. The air chooses at each position the easiest path between the grains, as a consequence of the competition between the air-water surface tension, the hydrostatic and the gas overpressure. In our case, we analyze the dependence on the two first quantities leaving the gas overpressure unmeasured.

We can separate the invasion process in two time scales, a short time scale and a long time scale, which we describe in the following.

Short time scale: Invasion regimes

At short time scale, we observe that the air invades the system in two qualitatively different regimes, which depend on the injection flow rate.

Low flow-rate: Percolation

• The air does not induce any significant displacement of the grains while creating different branches. The result of this process is the generation of a network that is more condensed near the inlet hole, while its density decreases closer to the free surface [Fig. 2.7(a)]. The emission at the surface is characterized by simultaneous branches. We can consider that, at least in the bulk, the granular material can be viewed as a simple porous medium2 .

High flow-rate: Fracture

• The air is able to displace the grains, thus fracturing the medium [Fig. 2.7(b)]. In a manner similar to the previous regime, the air creates several branches, but in this case, the dynamics is always governed by the higher branch. Consequently, there is one single emission point at the free surface.

As pointed out, the two regimes depends on the injection flow-rate. However, in an intermediate range of air flow-rate, both regimes where observed for the same set of parameters. Figure 2.7 shows the two types of air invasion, which correspond indeed to the same initial conditions. We can presume that the important variable behind this apparent behavior is the local rearrangement of grains (e.g. packing). As the system is polydisperse, and the packing can have slight differences from one experiment to another, the local conditions change, therefore inducing different regimes of air invasion.

Long time scale: Geometry of the fluidized zone

After the creation of the first air path the channel pinches and the air opens a new path. At long time, the air creates different successive paths and the channel network is permanently evolving. Due to the local rearrangements produced by these dynamics, the process leads to the generation of a fluidized region characterized by a lower compaction at the vertical of the inlet [Fig. 2.8(b)]. The fluidized region is centered at the injection inlet and widens while increasing in height z. The light transmission across the medium helps us to differentiate between the fluidized region and the exterior by producing a contrast between the lighter fluidized zone, produced by the lower grain compaction, and the darker region in the exterior. As in the previous 3D experiment, we aim at characterizing the bubble emission at the surface. But now, thanks to the formation of the fluidized zone we proceed with a slightly different analysis. We suppose (and it is experimentally corroborated) that the path followed by the air in the immersed granular layer is sensitive to local conditions only (e.g., grain polydispersity and packing), and not to what lies upon it (grain height). Therefore, it is equivalent to analyze the air passing through a virtual horizontal line inside the 2D granular bed (2D setup) to consider the bubble emission distribution at the surface of an immersed granular layer of total height z (3D setup). We thus measure the characteristic length ζ * (z) defined as the width of the fluidized zone for a height z. In order to determine ζ * (z), we let the system run for several hours, and measure the size of the fluidized zone inside the grain layer. This is done by image analysis, detecting the contrast between the two zones and then interpolating the results by a parabolic function3 [represented by the dashed lines in Fig. 2.8(b)]. We therefore get, in a single experiment, the characteristic length ζ * (z) of the region explored by the air as a function of the height z (Fig. 2.6). Similarly to the 3D experiment, we find that:

ζ * ∝ √ z , (2.4)
which is again independent of the flow-rate Φ. The result is compatible with the emission of bubbles at the free surface in a region having a typical size scaling like h g .

Venting dynamics

Comparison between 2D and 3D

From figure 2.6, we note that the slope of the two curves ζ(z) and ζ * (z) are different. This can be explained by two factors:

• First, as a consequence of the difference in the dimension (2D, 3D), we add in the 3D case one more degree of freedom (if we consider a square array). This gives the air more path to possibly explore, producing then, a slower increment in the height z which generates a lower diffusion coefficient.

• Second, we used two different methods to find the variable ζ and ζ * . On the one hand, ζ (3D) is taken as a typical size from a statistical analysis, whereas ζ * (2D) is defined as the total extent of the zone explored by the bubbles (maximum distance). We will not go farther on this difference, as the main purpose in this section is to point out the existence of the relationship stated in equation 2.3.

Analytical model

In this section, we develop a simple 2D analytical model to explain the empirical law ζ ∝ h g . The idea is not to explain quantitatively the experimental results, but to understand the basic mechanisms governing the dynamics of the system.

Description

We model the granular material by a rigid network of grains [Fig. 2.9(a)]. As mentioned in section 2.3 this could be the representation of the case where the system behaves like a rigid porous medium. The air path is generated when the inner overpressure (initially at the inlet) overcomes the capillary pressure associated with the void space between the surrounding grains. Due to the local heterogeneities (polydispersity or compaction for example), the path's upper end can propagate upwards or, more or less, sidewards. The results of this process is a meandering of the path towards the free surface.

We describe this process with a simple approach. We consider a square array where the air, at each step, can move to three different directions, left (-), right (+) and up (z) [Fig. 2.9(b)]. The propagation of air downwards is not allowed in the model. To take into account the local arrangement of grains we suppose that the capillary overpressure necessary for the air to overcome and propagate in the directions x < 0 (δP -), x > 0 (δP + ) and z > 0 (δP z ) can be written as:

             δP + = 2γ d δP -= 2γ d δP z = 2γ d -δP g (2.5)
where γ is the water surface tension and d the grains diameter. The three overpressures are not identical but taken from a random distribution around a typical value δP c .

The probability for the air to propagate upwards is biased by the hydrostatic pressure δP g ≃ ρgd/2 over the characteristic size of a grain d, where g is the gravitational acceleration, and ρ the water density.

The choice of this framework will neglect potential correlation between successive bubbles emission. Indeed as mentioned in section 2.2.2 this is experimentally corroborated, and in a first approximation, we can assume that there is no history in the system (the path that is generated will not bias the next path produced).

Probability for an upward propagation

Image analyses of the grains used in our experiments (in particular for d =250-425 µm) point out that the grains size distribution corresponds to a well-defined gaussian distribution (Sec. 1.4.1,Fig. 1.12). We can therefore fairly approximate the density probability for the overpressure by a gaussian distribution: .6) where δP is the capillarity overpressure around a typical value δP c and σ P the width which accounts for the local disorder in the material 4 .

Q(δP ) = 1 2πσ 2 P exp - (δP -δP c ) 2 2σ 2 P , ( 2
We want to find the probability for the air to choose the upward direction, denoted by p z . We use the notations to represent the probabilities p z (upward), p + (sideward,

48

Venting dynamics

δP g p (δP) σ [ p + , p -] [ p z ] δP 2 γ d 2 γ d δP g Figure 2
.10: Probability distribution for the system -The air probability for the three directions (p -, p + and p z ) are identical. They all share the same values for the mean µ and the standard deviation σ which account for the mean grain size and the associated error. Note that p z is biased by the hydrostatic pressure, which shifts the probability distribution, thus increasing the probability to find a lower value and hence increasing the probability for the air to propagate upwards. right) and p -(sideward, left). Figure 2.10 shows a graphical representation of these probabilities. We observe that p + , p -and p z have the same values of the mean µ and the standard deviation σ which represent the average grain size and the associated error respectively. As the system is subjected to the gravity, the upward probability p z is biased by the hydrostatic pressure δP g which shifts the probability distribution. The air will propagate upwards if p z is greater (δP z lower) than the other two equal probabilities p -(δP -) and p + (δP + ), we thus write:

p(δP z -δP g < δP + -) = 1 σ P √ 2π ∞ δPz-δPg exp δP -δP c σ P 2 d(δP ).
(2.7)

The three events are independent, and the conditional probability for the air to propagate upwards is governed by p z and [p(δp z > δp + --ρgd 0 )] 2 . Note that the second term is square, as it takes into account both the left and right directions. Finally, we can write the probability p z as:

p z = ∞ -∞ Q(δP z ) ∞ δPz-δPg Q(δP )d(δP ) 2 d(δP z ), (2.8)
which gives, for the gaussian distribution defined in Eq. 2.6:

p z (χ) = 1 π 3/2 ∞ -∞ e -v 2 ∞ v-χ -1 e -u 2 du 2 dv, (2.9)
where u and (v -χ -1 ) are defined by: (2.10) and the parameter χ is defined by: .11) The results show that p z (χ) is an integral function that it cannot be solved analytically. On the other hand, we have a function that depends on an unique dimensionless parameter (χ), which compares the width of the distribution of the capillary overpressure and the hydrostatic overpressure typical over a grain size. In the next section we build a diffusion coefficient that depends on p z and hence on χ.

u = δP + --δP c σ P , v -χ -1 = δP z -δP c σ P - δP g σ P ,
χ ≡ 2σ P ρgd . ( 2

Diffusion coefficient

We can estimate that the upper-end of the air path will move randomly N times to the left and right before moving upward, with N defined as:

N ∼ 1 p z -1.
(2.12)

As a consequence, the lateral mean-square displacement, associated with a displacement toward the free surface of the typical size d of a grain, can be estimated to: .13) We write the displacement of air in the x direction as ∆x = Ns(1 -p z ), where s is the unity of displacement (in particular, the size of grain d in the system) and N the number of steps. If we block the ẑ direction (p z = 0) all the steps N will be done in the x direction, i.e. ∆x = Ns; on the other hand, if we let the air pass only in the ẑ direction (p z = 1), we find ∆x = 0 . In the same manner, we can write the displacement in the ẑ direction as ∆z = Ns p z . Finally, we obtain a diffusion coefficient for a given upward probability (using d as the unity of displacement): .14) For a small (resp. large) grain size d 0 -→ 0 (resp. ≫ 1), the diffusion coefficient D -→ 1 (resp. 0). This behavior has been observed experimentally, and can be demonstrated numerically (section 2.5). Note that strictly speaking the coefficient that is defined here is not an usual diffusion coefficient. The units of an usual diffusion coefficient is [L 2 ]/[T ], which says, that the system evolves in time. As we pointed out before, we do not work with time and the system only evolves in space, which leads to a diffusion coefficient that depends only on [L].

δx 2 ∼ Nd 2 . ( 2
D(p z ) ∼ ∆x 2 ∆z = d 2 N(1 -p z ) dN p z = d 1 -p z p z . ( 2
p z = 1 -→ D = 0 p z = 0 -→ D = ∞ (2.15)
The results of the normalized diffusion coefficient D/d are shown in figure 2.11(a). We obtain a function that increases until an asymptotic value of D/d = 2 which corresponds to the maximum probability for going upward p z = 1/3. In the next section, we compare these results with the numerical case [Fig. 2.11(b), Sec. 2.5].

Numerical Analysis

We developed a numerical simulation with the same ingredients than the analytical model. In the simulations we can extract the data and compare the results with the previous case. This simple simulation is the starting point for future numerical simulations that will be extended in the next chapter.

Ingredients

Our first attempt to simulate the air paths consists in generating a path that depends only on the probability for the air to overcome the local overpressure. As in the analytical case, we choose randomly three values from a gaussian distribution, corresponding to the three different probabilities (p -, p + , p z ) that represents the capillary pressure. We control also the effective gravity in the simulation. This is done simply by subtracting a value (that represent the hydrostatic pressure δP g ) to the p z probability (see A global morphology is obtained, where a 'diffusion' profile z ∼ x 2 is represented by the red line. Note that even if the representation of p z = 1/3 means that there is no hydrostatic pressure, the choice of having three directions only generates an asymmetric invasion profile and not a circular (isotropic) profile as expected.

Fig. 2.10). Once we have generated the three probabilities we compare the values and choose the smaller, which represents the lower local overpressure and, hence, the direction in which the air propagates [Fig. 2.12(a)].

In the experiment, the generation of the fluidized zone is produced by the successive passing of air (after several hours). To mimic this in the simulation, and in order to get a good statistics, we generate 4000 different paths with an upper wall as a border condition. The results of the superposition for three and the 4000 paths are shown in figure 2.12(b-c). As we can observe, the superposition of all the paths produces in average a global morphology that is well-fitted by a parabolic function [represented by the red line in figure 2.12(c)].

Results

The most important quantities to compare in a diffusive system are the mean displacement ∆x and the square mean displacement ∆x 2 . To introduce these variables we use the Fick's law which describes the diffusion problem and can be used to calculate the diffusion coefficient D.

Fick's law

The experimental results suggest that we are in presence of a diffusion process. We can compare the numerical results with Fick's law which describes this type of system. Mean displacement for three different series. As expected all the data evolve around the injection point (dashed line, x 0 = 0). Note that the scale in the vertical axis is much smaller than the horizontal axis. (b) Evolution of the mean square displacement ∆x 2 for different values of χ. As in a random motion (diffusive system) the mean square displacement increases linearly in time (here in z). The arrow indicates the evolution of the parameter χ.

The flux density of the particles is proportional to the concentration gradient by: J = -D∇φ, (2.16) where J is the diffusion flux, D the diffusion coefficient and φ the concentration. Another way to represent the Fick's laws is the (second law): (2.17) which in 2D is analogous to the heat equation.

∂φ ∂t = D ∂ 2 φ ∂x 2 ,
We can use the density probability for the overpressure Q(δP ) as a solution for this equation. Resolving the equation we can extract two principal quantities that indicate if the system is diffusive or not: The mean displacement ∆x and the mean square displacement ∆x 2 . In a diffusive system, the mean displacement will be x 0 (the initial injection point). But the most important parameter will be the mean square displacement, which is a measure of the degree of fluctuations. For an initial concentration of the type Q(δP ), ∆x and ∆x 2 will be: ∆x = x 0 ∆x 2 = 2Dz.

(2.18)

From the numerical results we can calculate the mean displacement and the mean square displacement. These two quantities are represented in Figure 2.13. We can observe that the mean displacement fluctuates around the initial injection point (x 0 = 0) for every set of data. Figure 2.13(a) shows an example for three different values of χ. Note that the scale in both x and z directions are drastically different. We point out that, strictly speaking, a statistical average is taken for an infinite thermodynamic ensemble which is impossible in the experience. More important is the evolution of the mean square displacement [Fig. 2.13(b)]. As attended, similarly as Fick's law, the evolution of the quantity is linear in space (analog here to the time in a classical diffusion process).

If we take the slope of ∆x 2 vs z, we can estimate the diffusion coefficient D. The numerical evolution of D(χ) is displayed in figure 2.11(b) where we also compare the results with the analytical model. We find an excellent agreement between the data with only a difference in the scaling, maybe due to statistical reasons.

Discussion and conclusion

We have presented an experimental study of the gas emission location in 3D and 2D. The process is characterized by the propagation of air through an immersed granular bed and can be accounted by a diffusion-like model. We find the associated diffusion coefficient D/d = f (χ) which is a function of an unique dimensionless parameter:

χ = 2σ p ρgd . ( 2.19)
We predict that ζ * ∝ √ z which is compatible with the experimental observation ζ ∝ h g . The model predicts that the diffusion coefficient depends on the grain size d.

For small grain size (large χ), D increases, and the air explores a wide region around the emission hole (large D); for large grain size (small χ), the air propagates straight upwards and the region explored by the air is reduced (small D). Furthermore, σ p , which accounts for the distribution of the capillary overpressure, also depends on d.

Taking into account that the pore size (the space between the grains) is of the order of the grain size, one can estimate that the typical capillary overpressure δP c is of the order of γ/d, where γ stands for the surface free energy of the air-water interface.

Denoting by w the typical relative variation of the pore size, one can estimate further that σ p ∼ w γ/d and, thus:

χ ∼ 2wγ ρgd 2 (2.20)
Thus, we corroborate the dependence on the grain size of the control parameter χ and hence the coefficient D. For a granular bed made of monodisperse grains, the density can range, in practice, from that of the random loose packing, φ = 0.54, to that of the random close packing, φ = 0.64, so that w is of the order of a few percent. For a granular bed made of polydisperse grains, w is of the order of the relative width of the size distribution, about 25% in our typical experimental case in which d ranges from 250 to 425 µm. In this case, χ ∼ 10 and the model gives D ∼ 1.5d. Even if we cannot expect from such a simplistic model a quantitative agreement with the experiments, we can evaluate that, for a typical depth of h g ∼ 25 cm, the region of the free surface in which the bubbles are expected to exit must be of the order of dh g ∼ 1 cm, thus centimetric, in fairly good agreement with the experimental observations (Fig. 2.6). We remark that, for χ > 1, D/d remains of the order of the unity and ζ ∼ dh g .

Limitations of the numerical analysis

As presented in the chapter, the numerical simulations and the analytical model are based on strong hypotheses:

• The numerical simulation (as well the analytical case) does not allow the downward propagation and even do not compare in each moment the neighbors probability. It is evident that this supposition is not realistic. Also, we only calculated the probability in one point (the air path tip) and not between all the neighbors (as in a the real experience), neutralizing the generation of branches. This will be included in the study in the next chapter. In the next section, we add more ingredients to the simulation with in particular, the possibility of branch formation.

• The dynamics is not taken into account neither in the analytical model nor in the numerical simulation. For example, assuming a rigid network can never lead to the fracture regime observed in the experiments.

• As we chose at each time a new probability, the local history of the path was not considered. It means that from time to time, the air can occupy the previous position which does not conserve the total volume of the injection.

• The history was also not taken into account. When a new path is created, the next path is created from scratch. As we know from the experiment, the path that is created changes locally the medium properties, and thus the dynamics of the next path. The latter could then exhibit a different dynamics that we fail to describe with this simple model.

Cluster formation

An interesting phenomenon was observed from time to time in the 3D experiment.

Typically for high values of h g (h g > 6 cm), we observe that the bubble emission remains in a very localized region forming clusters of bubbles. The results for the bubble position r(t), show occasionally the stabilization of r around a constant value for several bubble emission revealing the formation of such cluster (Fig. 2.14). In this case, the air propagating through the immersed granular layer is able to create and maintain an open channel, over most of the height, and this path remains open for several bubbles passing -with small variations of the bubble emission location at the surface. This type of regime (open channel) has been already observed in a similar system [START_REF] Gostiaux | Dynamics of a gas bubble rising through a thin immersed layer of granular material: an experimental study[END_REF]. Here, we point the difference with previous observations, i.e., the existence of such channel in the granular layer, leading to cluster formation at the surface. This phenomenon might be due to local fluidization, which makes possible the small variations in position. Note that no cluster is observed for small values of h g . Indeed, the air explores a region too close to the center, and the only possible 'cluster' is the centered distribution of the bubbles around the vertical of the injection point.

One key indicator for this phenomenon could be the autocorrelation function mentioned in section 2. 2 [Fig. 2.5(b)]. Indeed, in the previous study, no long term correlation was observed in the signals. But the dependence of this parameter as a function of the problems parameters (i.e. flow rate Φ or water height h g ) was not studied. Finally, as the methods of observation where different in the two setups, we cannot quantify the apparition of cluster in the 2D system. Nevertheless, the branching formation in both regimes (fracture, percolation), and the generation of the fluidized zone could be important indicators of the possible appearance of clusters.

Introduction

In the previous chapter, we showed that the air invasion through an immersed granular bed can be interpreted in terms of a diffusive-like model. We proposed that the diffusion coefficient D is an increasing function of a unique parameter of the system, χ ≡ 2σ P /ρgd which compares the width σ P of the distribution of the capillary overpressure associated with the passage between the grains with the typical hydrostatic pressure variation ρgd over the grain size d. The next logic step would be to observe in the experiment the dependence of D on χ. Moreover, it would be necessary to have a more realistic theoretical approach. If we take a closer look at the parameters governing χ = χ(σ, ρ, g, d) we observe that from a practical point of view almost all of them are not fit to be used. The capillary pressure variations σ, the fluid density ρ and the grain size d demand at each time a new experiment which is time-consuming1 . However, the only parameter left, the gravity g, can be controlled by tilting the experimental cell, which is easier and faster. This allows us, without changing the medium, to span almost all the range of χ, from a vertical cell where g eff = g (χ = finite) to an horizontal cell where g eff = 0 (χ = ∞).

Motivated by the observation in nature such as the gas emission at the sea floor, the oil leakage from active hydrocarbon systems (pockmarks) or a deep overpressurized petroleum reservoir (Sec. 1.1.1) we extend the previous numerical work. We then will concentrate in the study of a constant volume of air injected and the potential spatial extension which could be helpful in determining the potentially invasion zone.

In this chapter, we analyze the morphology of the invasion region generated by successively injecting a constant air volume observing the formation of a region with a delimited contour. We contrast preliminary experimental results with numerical simulations. The simulations are divided in three steps (Sec. 3.3.1). First, we observe the geometrical characteristics of the invaded region. We then observe the evolution of the density profiles and use these results to observe the evolution of the geometry as a function of the dimensionless parameter χ and the volume injected. Finally we define the aspect ratio of the invasion zone (Sec. 3.3.2). We then change the bottom boundary condition and compare the results for both simulations (Sec. 3.4). We analyze the extension and velocity of the invasion zone (Sec. 3.4.1) and also its curvature closed to the injection point (Sec. 3.4.2). Finally we define an effective diffusion coefficient (Sec. 3.5).

Preliminary experimental results

We performed a series of experiments in the 2D Hele-Shaw cell (see Sec. 2.3), observing the air invasion pattern when injecting a constant volume V . The results for a fixed volume of air injected show that for the vertical configuration (g eff = g) the air invasion is an almost straightforward path without much branches [Fig. 3.1(a)]. Tilting the experimental cell, we reduce the effective gravity, and for a completely horizontal configuration (g eff = 0) the branches are more numerous and, thus, the pattern is more compact [Fig. 3.1(b)]. Accordingly, the invaded region is broader and less high in this case whereas an increase of g eff leads to a narrower and taller pattern with less branches. In order to characterize the morphology of the region potentially explored by the air for a given volume injected, we repeat the experiment, and superimposed the different air paths obtained for the same initial condition2 . In order to get a good statistics and to determined precisely the contour of the invaded region, we would have to repeat the experiment a large number of times, which is time-consuming. Indeed, once the pattern is formed, in order to reset the initial condition, one must open the cell to remove the air trapped in the system, which takes some time (Sec. 2.1.2). Moreover, we would not be sure to prepare the system in the exact same way and the state of the granular packing (its density, for instance) might be different from one run to another. As a consequence, we concentrate in the numerical work adding more features that approaches to the real experiment.

Model description

Ingredients

The previous numerical model (Sec. 2.5) only took into account the probability to propagate in three directions (up, left, right) at each point, while the downward propagation was prohibited. Furthermore, as a consequence of choosing at each time a new probability, there was no local history of the path (the air can circulate over two adjacent points without extending the path). Here, we improve the simulation by changing the basic ingredients of the model.

To take into account the local history of the path, in a more realistic model, we simulate the experiment by generating a bidimensional square network where each vertex of the cell represents the void space between four grains. The size of the network is about 1200 × 1000 grains, that is nearly the same in the experiment (around 1600 × 1200) 3 . This is sufficient for the air to avoid the lateral borders without affecting the system (as we are interested in quantifying the morphology of the air invasion we do not allow the air to reach the surface.). The only condition that we impose is the restriction at the bottom wall, through which the air cannot propagate. This can affect the air path, specially when χ ≫ 1.

Initially, each value of the network is completed with a random value obtained from a normal distribution [Fig. 3.2(a)]. This value (with a standard deviation of σ p = 1) is related to the grain size and hence to the capillary pressure given by σ P = 2γ/d where γ is the air-water surface tension and d the grains diameter. Once our network is filled, we proceed by adding a vector that represents the hydrostatic pressure δP g = ρgd/2, where ρ is the density of the water, g the acceleration due to gravity, and d the typical grains size, as a control parameter as in the tilted experiment. When the cell is completely horizontal, there is no hydrostatic pressure (δP g = 0) and when the cell is completely vertical δP g = ρgd/2. Finally, we consider the dimensionless parameter that governs the dynamics of the system χ = σ P /δP g = 2σ P /ρgd, ranging from 10 0 (vertical cell) to 10 4 (horizontal cell).

Protocol

We start from the first point at the bottom center of the network (corresponding to the virtual air injection point) and compare the local pressure values of the three first neighbors (left, up, right). We then choose the lower value which corresponds to the easiest path for the air to propagate. Because the initial position is at the bottom of the cell, at this step the air cannot search for a path below its position. We repeat the process by comparing for each step all the neighbors for every position along the path and choosing the lowest one. With this type of condition, the model allows the generation of branches [Fig. 3.2(b-c)], process also observed in the experiment (Fig. 3.1). Because the idea is to observe the path of air for a limited volume, we limit the exploration of the system to a fixed number of steps (N = 2000, 4000 and 8000 in the simulation) as an analogy for a volume fixed experiment. The volume can be defined as V 0 = Ns p d where s p is the surface associated with a pore and d the grains diameter. The final invasion geometry is obtained by the superposition of 2000 simulations for each set of parameters (χ, N).

Aditional simulation: Bottom boundary condition

In order to reproduce the experimental condition, the first numerical simulations were performed with a bottom wall, prohibiting the air propagation below the altitude of the injection point. However, in order to account for the difference with the air injection in a free medium, we also performed simulations without a bottom wall, thus allowing the air propagation below the injection point (Sec. 3.4). The ingredients and protocol are the same in both simulations, using the same set of parameters for χ and adding one more volume N = 6000 for this last case.

Morphology of air invasion: Numerical analysis

Morphology of air invasion: Bottom wall

In the following we describe the geometrical characteristics of the invaded region and its dependence on the volume N and the dimensionless parameter χ. In order to quantify the morphology we characterize the contour of the invasion obtained by the density profiles. Finally, we use the typical width and height to characterize the invasion region, and investigate the evolution of the ratio between these two quantities (aspect ratio) as a function of χ.

Geometrical characteristics of the invaded region

The superposition of 2000 numerical paths makes it possible to observe the geometry of the global invasion pattern (Fig. 3.3). We observe a drastically different geometry that depends on the dimensionless parameter χ that we describe in the following.

Qualitative observations

• As expected for χ = ∞ [δP g = 0], the air invasion is almost isotropic [Fig. 3.3(a)].

The air has no preference for a particular direction and in average travels the same distance in every direction, thus forming roughly half a circle.

• Adding a pressure gradient (δP g ) to the system, we break the symmetry of the air invasion. Decreasing χ (increasing the gravity) the invaded region elongates along the vertical. The resulting geometry is close to a parabolic shape that has almost the same width while the height increases when lowering χ [Fig. 3.3 (b-e)].

Nevertheless, the simulation presents an anisotropy due to the presence of the solid bottom wall. To take this into account, we analyze the evolution of the density, which gives another manner to quantify the morphology.

Density and probabilities profiles

As mentioned before the invasion is not homogeneous in space producing a symmetry breaking in the density. This is caused by two factors. On the one hand, as a consequence of the bottom wall, the exploration is denser near the bottom, increasing the number of time that air crosses a vertex, specially for χ ≫ 1 (lower g eff ). On the other hand, the injection is realized at one single point, and thus, the density is the highest around the injection point. Figure 3.4 shows that invasion density for χ = 125, and N = 8000. In Fig. 3.4(a) we observe a surface plot for the normalized density ρ p = ρ/2000, where ρ is the number of time of a section has been crossed. At the initial injection point (x i , z i ), the normalized density reaches its maximum value ρ p = 1 and as we move around (more or less 20 vertex) in any direction, ρ p decreases to a value of ρ p = 0. 25 [Fig. 3.4(b)]. In summary, we observe that almost all the invasion takes place for a density lower than ρ p < 0.25. We can define the function f iso which represents the normalized percentage of number that the section has been crossed. For example f iso = 1 % means that the vertex has been crossed 2000 × 0.01 = 20 times.

From the invasion densities we can obtain the profiles at a given height z [Fig 3.5]. This quantity is helpful in determining the typical width of the invasion region. Figure 3.5(a) shows three normalized density profiles at different heights z. We observe that the distributions go from a delta type when z = 0 (injection point) to a normal type distribution for higher values. We fit the results with a normal distribution centered at the initial condition x = 0 and take the values of the standard deviation σ that account for the width of the distribution and hence the width of the air invasion at a fixed height z. From the profiles4 , one gets the typical width σ of the region crossed by the air paths as a function of z [Fig. 3.5(a)], which gives us the evolution of the mean square displacement as a function of height [Fig. 3.5(a)]. We observe that for small enough values of χ and z, σ 2 is proportional to z, which accounts for the parabolic shape of the invaded region near the intlet.

Morphology of air invasion: Numerical analysis

We also analyze the evolution of the invasion geometry qualitatively. We determine the contour of the invasion zone by finding the isodensity lines. In order to do so, we search from the sides toward the center a fixed value that represents the number of times that air has passed in our superposition path's matrix [blue lines in Fig. 3.7(a)]. This makes it possible to observe the evolution of the invaded zone for the different parameters.

Dependence on the injected volume N and χ

Observing the isodensity lines for f iso = 2 %, we note that the geometry of the invasion zone is highly dependent on N and χ. From Fig. 3.6(a), we observe the contour plot for a same χ and for four different volumes injected. Contrary to the intuition, the growth of the invasion zone is not a simple dilation, instead, we have a faster growth in the vertical direction. This observation is more evident looking at the evolution, for a same volume injected and different χ [Fig. 3.6(b)]. As χ decreases (increasing δP g ) we transit from a non-diffusive system (isotropic invasion) to an almost diffusive behavior (parabolic contour). The limit is reached for χ > 2, in this case, the air goes straight up forming a vertical line.

Aspect ratio of the invasion zone

In order to quantify the evolution of the invasion geometry as a function of N and χ we define the typical width and height on both axes ∆x and ∆z. To do so, we consider the intensity profile n x (x) (reps. n z (z)) along the x (reps. z) axis. The intensity profiles are defined by the addition of all the number of times that air has passed over the whole height [n x (x), Fig. 3.7(b)] and width [n z (z), Fig. 3.7(c)]. As expected, n x (x) is symmetric around the injection point, and we consider ∆x to be the width at half the peak value. On the other hand, n z (z) [Fig. 3.7(c)] is almost flat over the whole height and suddenly decreases above a given altitude. We define ∆z as the altitude of the point corresponding to half the plateau value. Now that we have the two typical width and height, we quantify the invasion geometry using the aspect ratio ̺ = ∆z/∆x as a function of χ [Fig. 3.7(d)]. The curve is globally the same for all the data for different injected volumes. We observe that ̺ ∝ 1/ √ χ in a wide range of χ and reaches a plateau value of the order of unity (̺ ≃ 0.82) for large χ. For a compact structure (isodensity), one would expect the aspect ratio to tend to 1 in absence of gravity, in absence of bottom wall. Taking into account that the density is not constant within the structure and that the air flow is limited by the lower edge, we are not surprised that the ratio tends to a value of the order of unity but not exactly5 to 1. Interestingly, we note a collapse of the data when reporting ̺ as a function of χ/ √ N : the scaling ̺ ∝ 1/ √ χ holds true for χ 30 √ N whereas ̺ ≃ 0.82 for larger values.

Air invasion in a free medium: Comparison

In order to investigate the effect of the bottom wall which limits the invasion, we performed the simulation in a free medium (i.e. without a solid bottom wall). In the following, we compare the results for both boundary conditions (with and without a bottom wall) presenting numerical results for the invasion in a free medium for 2000 simulations for the same set of parameter χ and injection volume N = 2000, 4000, 6000, 8000. We contrast the values for the typical height ∆z and width ∆x, and compare the normalized aspect ratio for both boundary conditions. We calculate the invasion velocity in both directions as a function of the injected volume. Finally, we quantify the curvature of the invasion zone close to the inlet, which gives us the possibility to explore the type of function behind the invasion process. First, we compare the values of the typical length ∆z and width ∆x for the two bottom conditions (Fig. 3.8). As expected we find that the geometry of the invasion zone for χ = ∞ is well defined by a circle with the ratio6 ̺ = 1.03. In contrast with the same χ but with the bottom wall condition we find that the ratio is ̺ = 0.82 [Fig. 3.8(up)]. Comparing the two values for ∆z and ∆x we find a difference of the order of 10% between both condition for χ = ∞ and a difference of 2% for χ = 250 [Fig. 3.8(down)]. The difference is due to the fact that for high χ the bottom wall plays an important role restricting the air path. When adding the pressure gradient δP g the vertical axis is more predominant, leaving the bottom unaltered for both boundary conditions.

We present the aspect ratio ̺ = ∆z/∆x vs χ/ √ N for the two conditions (Fig. 3.9). From the figure we observe that both data collapse over almost all the range in χ/ √ N having the same slope: ̺ ∼ χ -0.59 which is almost proportional to 1/ √ χ. We find a slight difference between the two conditions for χ/ √ N 30.

68

Morphology of air invasion: Numerical analysis

Extension and velocity of the invasion pattern

As reported in section 3.3, it was qualitatively observed that the variations of the typical height ∆z and width ∆x were not similar. Globally, we observed a linear increases in ∆z and a slow increases in ∆x as a function of the injected volume. Also, as the pressure gradient in the system goes higher (lower χ) the invasion is faster in the vertical direction. Reporting the values of ∆z and ∆x as a function of the volume [Fig. 3.10(a-b)], we observe that the evolution of the displacement ∆z [Fig. 3.10(a)] is linear in space, its slope increasing while decreasing χ (higher the pressure gradient, the faster the growing). On the contrary, the displacement ∆x [Fig. 3.10(b)] is an unknown function that has a linear regime that gets slower when increasing the volume. The growing of ∆x increases, while increasing χ (the smaller the pressure gradient, the faster the growing).

We can extend the results and calculate the normalized velocities, taking the slope at the origin, in both directions (see insets Fig. 3.10). It is interesting to notice that the horizontal speed ν x /d is almost constant in χ around a value of ν x /d ≃ 10 -2 . In contrast, ν z /d is a decreasing function of χ reaching the same asymptotic speed ν x for high χ. This is the case when there is no gravity, and the invasion is isotropic (no wall): as expected, both velocities are the same. We can find a two orders of magnitude difference for lower values of χ, which explains the great difference in the aspect ratio of the invasion.

It is clear that the displacement in ∆z is linear in volume, but it is not the case for ∆x. In order to investigate more precisely this behavior, we analyzed in more details the evolution of ∆x and ∆z as a function of N, for a fixed value of χ (Fig. 3.11). As expected, we find again the linear evolution for ∆z [Fig. 3.11(b)]. The results for ∆x

Curvature of the invasion zone

The invasion profiles without a bottom wall make it possible to compare for different sets of parameters the curvature of the invaded region close to the inlet. The simulation with the bottom condition cannot let us have access to this curvature, and we can only find that the invasion region has a quadratic profile z ∼ x 2 (Sec. 3.3) indicating that the curvature, at least at the beginning, corresponds to a parabolic shape (diffusive process).

To investigate the initial radius of curvature, and hence, determine if the invasion process is diffusive or not, we define the distance l f that separate the injection point (x i , z i ) with the base of the resulting geometry obtained for the isodensity lines f iso = 1 % [represented in Fig. 3.12(a)]. The parameter l f is defined as the distance that separates the injection point (x i , z i ) with the bottom base of the resulting geometry. In red is represented the initial contour for f iso = 1% which is fitted by

f (x) = x 2 /2a. (b) Ratio of l f /R c (0), with R c (0) = a (see text).
The results show a transition at χ ∼ 10 3 . Over this value the invasion is close to a circular (non diffusive) invasion while under this value the parabolical shape accounts for the diffusive behavior.

Diffusion coefficient 71 Initial Radius of curvature

The radius of curvature of the invasion profile can be the critical parameter to differentiate if the invasion is a diffusion process or not. It is defined by:

R c (x) = (1 + f ′ (x) 2 ) 3/2 f ′′ (x) , ( 3.1)
where f (x) is the function fitting the invasion region contour and f ′ (x) and f ′′ (x) its first and second derivatives, respectively. Supposing the simple case of a parabolic function centered at the origin defined by f (x) = x 2 /2a and calculating the first and second derivatives, we find that:

R c (x) = 1 + x 2 a 2 3/2 a . (3.2)
As we search for the initial curvature, we can simplify the problem by considering only the curvature at the center of the parabola. This gives us:

R c (0) = a (3.3)
which means that the radius of curvature at the origin is directly given by the coefficient of the parabola. We extract the value of a, by fitting f (x) = x 2 /2a to the invasion contour with a value of f iso = 1% for the isodensity lines [Fig. 3.12(a)].

Comparing the ratio of the characteristic length l f with the radius of curvature at the origin R c (0) we observe two well-defined regions for the diffusive and non-diffusive regimes [Fig. 3.12(b)]. A transition point at about χ ∼ 10 3 is found. When the ratio goes to 1 it indicates that both lengths are equal and then, the invasion is isotropic. Surprisingly, for lower values of χ the ratio does not goes to 0 but forms a plateau around l f /R c (0) ∼ 0.3.

Diffusion coefficient

We now know that the system is diffusive for χ 10 3 . It is appropriate then to find the associated diffusion coefficient. We process by two methods, the first, accordingly to Fick's law, by using the evolution of the mean square displacement ∆x 2 vs z, and the second, by characterizing the ideal case of an air path that crosses the system.

Mean square displacement ∆x 2

In Chap. 2, Sec. 2.5 we found the diffusion coefficient by calculating the mean square displacement ∆x 2 for the numerical simulation. We can proceed in the same manner with the results of the superposition of the invasion paths (Fig. 3.3). From the intensity profiles we have already calculated σ 2 [Fig. 3.5(b)]. According to Fick's law, taking the slope at the origin gives us the associated diffusion coefficient 7 . The results of D/d vs χ are shown in figure 3.13(a). We found that D/d scales like ∼ χ 2/3 and as expected there is a transition around χ ∼ 10 3 , above which the invasion is not diffusive anymore.

Effective diffusion coefficient

If we considerer the regime where the system is diffusive (low χ and large volume injected), we can define an effective diffusive coefficient. We consider the ideal case of a path traveling in the system with a parabolic movement z ∼ x 2 . Then, the total length of the arch is given by ds

= 1 + f ′ (x) 2 . Using f (x) = x 2 /D eff (diffusion-type process), we have: l = x 0 0 1 + 4x 2 D 2 eff dx , (3.4)
where x 0 is the position of the upmost point of the path on the horizontal axis. We can calculate the integral considering the case for large volume injected x 0 ≫ 1:

l = x 2 0 D eff . (3.5)
For x ≫ 1 the arch length l is the same as the length of a straight path going upward. We can then define the aspect ratio ̺ for a high volume injected as the ratio of the two principal lengths:

̺ = l x 0 = l D eff . (3.6)
The normalization of the aspect ratio can be obtained by considering the injected volume V 0 = s p l, where s p is the unity pore surface. We then obtain:

̺ √ V 0 ∼ ̺ s p • l = 1 s p • D eff (3.7)
From this last equation we can define the effective diffusive coefficient given by D eff ≃ d χ provided that the gravity effects are large enough (χ 30 √ N):

D eff = V 0 s p 1 ̺ 2 (3.8)
Using V 0 = Ns p d, we obtain an effective diffusion coefficient as a function of the volume (numbers of steps N) and the aspect ratio ̺:

D eff d = N 1 ̺ 2 (3.9)
The evolution of the effective diffusion coefficient as a function of χ is reported in Fig. 3.13 (b). The collapse of all the data is obtained by normalizing the aspect ratio ̺ by the injected volume N.

Discussion

The model presented here has a richer behavior than the exhibits in the previous analytical and numerical models (Chap. 2). In the previous model, branching formation and the downward propagation of air were not allowed. This last condition is fundamental to understand if the system is diffusive or not. Indeed, a diffusive system is not reversible in time. Here in the simulation, time is analog to the height z. As the system can generate branches and also search for a path for lower z, the system does not satisfy the diffusion condition. The condition is only satisfied for lower values of χ, or in other words, for a high hydrostatic pressure δP g for which air can almost only propagate in the vertical direction. We showed however, that even for a non-diffusive behavior, we can define an effective diffusion coefficient.

Relationship between χ and branches formation

It is particularly interesting to interpret the parameter χ ≡ 2σ P /ρgd, which compares the width of the distribution of the capillary overpressure within the pores to the Morphology of air invasion: Numerical analysis variation of the hydrostatic pressure over the grain size. On the one hand, for a given injected volume N, the air inside the existing path being connected, one can consider the path as an isobar. On the other hand, the pressure in the liquid, outside the air path, increases from the path tip (its highest point) downwards, which makes the formation of a side branches more and more difficult at depth, below the path tip. Considering the physical meaning of χ, one can estimate that side branches cannot form at a distance larger than d χ below the tip. As a consequence, the gas injection results either in the formation of side-branches in a region of typical height d χ (in a local increase of the gas fraction or in the widening of the structure) or in the growth of the tip upwards. The result of such complex dynamics is a complex path exhibiting more or less side branches depending on the value of χ. Interestingly, in average, the gas occupies an elongated region whose radius of curvature near the inlet is d χ, as shown by the numerical results.

Gaseous fraction and maximum height

From the dependence of the height ∆z and typical width √ D eff ∆z, we can estimate the typical gaseous fraction F inside the invaded region. Estimating the corresponding surface area √ D eff h 3/2 , one gets from simple algebra

F ∼ χ 1/4 / √ N , (3.10)
thus dependent on χ and on the injected volume (on N). The fraction F slightly increases when χ increases, i.e. when the gravity effects are reduced and the side branching enhanced. In addition, denoting w the typical relative variation of the pore size as proposed in the previous chapter, one can estimate further that σ p ∼ wγ/d and, thus, that χ ∼ wγ/(ρg eff d 2 ). Thus, for a given (natural) system, considering that the pore volume v p ∝ d 2 and taking into account the result obtained for ν z , we are taught that the maximum height z max reached by the gas within the granular does not depend on the grain size and scales like:

z max ∼ V /(l c √ w) (3.11)
where l c ≡ γ/ρg eff . Thus, provided that the proposed estimate of χ is correct, for a given volume V , the maximum height is controlled by the capillary length l c and the relative width w which accounts for the heterogeneity of the capillary overpressure.

Conclusion

In conclusion, we reported results of a numerical study which makes it possible to predict from the knowledge of one single control parameter χ, the typical height, width and gaseous fraction of the region invaded by a given volume of gas liberated at the base of an immersed granular bed. Such results could be of practical importance: for instance, gas can be trapped on purpose in an underground natural container below a sediment layer. Our results can help in predicting if the gas is likely to reach the free surface and escape the system if the reservoir presents a defect (hole or fracture).

We propose in a future work, to extend from the theoretical point of view, to slightly different geometrical situation, especially to the 3D case and, from the experimental point of view, to the case of an horizontal liquid flow.

In the previous chapters we have focused on the final state of the global invasion pattern in 2D and the zone explored by the air in a 3D configuration. The numerical results obtained indicate a linear growth in z and square root in x. In both cases the model used to describe the evolution of the system was a static model, where the history and the inertia of the air flow were not taken into account.

In this section we present preliminary experimental results of the evolution of the two interfaces delimiting the fluidized zone: (1) in the bulk, the interface between the grains in the fluidized zone and the motionless grains, and ( 2) at the surface of the granular bed, the interface between the fluidized zone and the water. First, we focus on the geometry at long time of the fluidized zone for different flow-rates (Sec. 4.1.1) and quantify the evolution of the diffusion coefficient as a function of this parameter. We then repeat the process for different angles θ by tilting the cell (Sec. 4.1.2) which is similar to vary the dimensionless parameter χ. After describing the stationary state, we account qualitatively for the dynamics of the bulk interface and the mechanism involved (Sec. 4.2.1). We estimate, in particular, the speed of the interface front for a small region of the space (Sec. 4.2.2). Finally, we describe the evolution of the free surface of the granular bed (Sec. 4.3).

Stationary state: Long time scale dynamics

As already presented in section 2.3.2 the fluidized zone is generated by the successive paths of air through the granular bed. Nevertheless, the effect of the air-flow in the system was only studied for the initial injection regime (Sec. 2.3.1) and not for a longtime experience. We might think that changing the flow-rate can affect the final form of the fluidized zone.

In order to quantify experimentally the dependence of the fluidized zone morphology on the flow-rate Φ and the angle θ, we analyze the stationary final state for the injection of a large volume of air. Contrary to the previous numerical case (Sec. 3.3), where the air-path was sufficiently small not to reach the surface (as we were interested in the morphology of the invasion zone), here, we inject a sufficiently large volume of air to form the fluidized zone (we inject air at a constant flow-rate Φ for typically a few hours, up to one day). For the different flow-rates, the total duration of the experiment is chosen so that the same volume of air has crossed the system.

Dependence on the air flow rate Φ

The final state is obtained by adding all the images for one series for a fixed flow-rate. This method provides the best results for the latter analysis 1 [Fig. 4.1 (a)]. From this image we detect the contour of the fluidized zone (black lines). In order to avoid the slightly disperse around an average diffusion coefficient D = 2.12 ± 0.15 √ cm.

Based on this preliminary results, we can conclude that, al least for θ = 0 • and the range of Φ used, the history and the inertia do not matter in the generation of the fluidized zone. We do not find any significant variation of the diffusion coefficient D as a function of the flow-rate Φ for a same total injected volume.

Dependence on the tilting angle θ

In a similar way, we repeat the process for different angles θ and find the diffusion coefficient. The results shows a non trivial behavior (Fig. 4.2). First, we note that we can only analyze the data for values lower than θ ≃ 70 • . Above this angle, the channel, opened by the air, stabilizes due to the fact that the air pushes the grains out of the initial granular bed. The gravity is not able to close the channel, which is therefore stable in time. Note that in this case, the grains do not occupy anymore the whole gap of the cell. This is dramatically observed for values of θ higher than 70 • [Fig. 4.3]. From the images we show two examples for θ = 80 • where we find a straight air path [Fig. 4.3 (a)] and a zigzag path [Fig. 4.3 (b)]. Both images are taken right after the injection starts. After some time, the system expels the grains, forming a delta (like rivers), and the only material left in the center of the air path is the water.

The dependence of the diffusion coefficient D on the angle θ is not simple to interpret. In the previous section, we found that the effective diffusion coefficient was an increasing function of the parameter χ and hence a decreasing function of the gravity g. The results shows a completely different behavior. Here, we have a slightly decreasing evolution of D as a function of χ(θ) with an apparent abrupt change at about θ ∼ 45 • . The assumption of the dependence in the injected volume only leaving any effect of the inertia can actually change the dynamics, specially when the cell is tilted and the inertia is predominant. Maybe one possible explanation is that the problem depends only on its geometry (grain size, gap space, inlet diameter, etc) and not on its dynamics. A similar problem has been study by [START_REF] Melo | On drawbody shapes: From Bergmark-Roos to kinematic models[END_REF]. In that case, the problem is inverse, they studied the geometry formed by the ore recovery in underground mining called draw body shape. The results indicates that the final geometry depends mainly on the volume of the extracted material.

Transient state: Bulk interface

The transition to the final state of the fluidized zone was only studied numerically, observing its evolution for a number of fixed steps (Sec. 3.4.1). The dynamics inside this region as well as its growth and speed were not studied experimentally. In this section we present preliminary observations and we describe qualitatively the mechanism behind this formation. We also analyze the velocity of the interface motion during the growth of the fluidized zone.

Bulk deformation

Here, we present preliminary observation of the transient dynamics of the interface separating the fluidized zone from the outer region where the grains remains still.

• When the air starts crossing the system, a fluidized zone is generated first at the free surface of the granular layer, deforming locally a region closed to this boundary. Then, in time, this region propagates toward the bottom [Fig. 4.4 (up)]. This can be explained by the fact that the system is not constrained in all the directions and the loose region close to the surface is moved more easily.

The borders of the fluidized zone are well-defined due to the change of density (compaction) between the interior (fluidized zone) and the exterior (compact zone).

• Inside the fluidized zone, we report an interesting dynamics which can be observed thanks to the bubbles which are trapped inside the bulk. Indeed, they act like tracers indicating the direction and speed of the grains flow [Fig. 4.4 (down)].

First, the grains are pushed away from the incoming air path. As this process evolves constantly in time, some bubbles get trapped at the interior (and exterior) of the fluidized zone. Surprisingly, through time, these bubbles do not rise up to the surface but descend toward the injection point. The speed of the descending bubbles depend on their horizontal position, getting slower when close to the interface. In the process, bubbles can be created (by the air path), coalesce with others bubbles or get caught by the ascending air path. This type of movement can be related to a convective process related to the grains flows.

t = 0 s t = 1h 15min t = 2hr 45min t = 18hr 15min 5 cm t i Δt = 1hr 2 cm t f
From this latter observation we must point out that we are really in presence of a 'fluid' region that evolves permanently in time with specific characteristics.

Interface speed

From the numerical simulations we have obtained the growth and the velocity of the global invasion pattern (Sec. 3.4.1). In that case, the results were obtain by taking the characteristic lengths (∆x, ∆z) along both axis. To obtain the experimental global profile at each time is complicated. For one part, the growth is not homogenous in both directions (the fluidized zone starts at the free surface and then stabilizes) and for the other part, we do not have the spatial resolution. We then choose a different method. We concentrate on a specific small region near the center of the cell and analyze the evolution in that zone of the fluidized front.

In order to quantify the growth of the fluidized zone, we follow at different positions (horizontal and vertical) the fluidized interface [Fig. 4.5 (a)]. This latter is defined by a well-marked contrast (darker grains) 3 . The results for the horizontal position front L f (x) vs t are shown in figure 4.5 (b). We find a fast growth at short time scale which gets slower in time. The data are well-fitted by a logarithmic function of time for the different heights. When normalizing the data by the square root of the height z (we must remember that the borders of the fluidized zone corresponds indeed to x vs √ z) we observe a collapse of all the curves [Fig. 4.5 (b) inset]. The movement of the interface front is not homogenous in time, for the first minutes we observe a strong acceleration, and then its stabilizes. Additionally we report the evolution of the vertical velocity [Fig. 4.5 (c)]. In this case, the velocity of the front is faster in the vertical direction than in the horizontal direction. These observations, and the presence of small events in the evolution of L f (x, z) indicates that there are two mechanisms of grains transport.

We notice that close to the fluidized region, there is a well-defined zone which displays spontaneous events in time. These events are similar to the stick-slip dynamics [Fig. 4.6]. This stick-slip zone (SSZ) evolves in time, and is linked to the fluidized zone. In fact, in time, the fluidized zone increases slowly while the SSZ gets thinner and finally disappears after some time. We can understand this type of dynamics as a consequence of the air path, which is continuously shearing the bulk producing a deformation that (we suppose) generates the stick-slip motion. The combination of this dynamics plus the effects of the grains flow can possibly by at the origin of the faster growing at early stage and the small events reported in L f (x, z). In any case, the modeling of this type of movement is not easy and we restrict the analysis only to this.

L f (x) (cm) L f (z) (cm) L f (x) (cm) (b) (c) (a) √z x 1 x 3 x 2 z 1 z 3 z 2

Dynamics of the free surface

In time we observe that the evolution of the free surface is linked to the fluidized zone dynamics. It has been shown (Sec. 4.2.1) that the formation of the fluidized region starts indeed at the free surface and then propagates toward the bottom. Here, we analyze the free surface motion in time and quantify the emission of bubbles along this surface, above the fluidized zone. The free surface is detected by image analysis. We considerer three different quantities (Fig. From the figures we observe that the dynamics is concentrated around the center of the cell. As expected the borders remains constant around the initial injection point through time (except from some isolated points due to the initial 'percolation regime' invasion). The variations of height [Fig. 4.7 (b)] can help to locate the exact position of the emission showing that sometimes the air-path stays in a fixed position 4 . Binarizing this last result [Fig. 4.7 (c)], makes clear that the exploration over the free surface does not change in time in the stationary regime (the red color represents bubbles emission while the blue color indicate static zones). Adding all the numbers of events for each column and row, gives us the numbers of bubble emitted at each position and time, respectively. We report a centered distribution around the inlet with a flat peak5 . Note also, that the emission of bubbles is constant in time.

This may suggest that the maximum range of the air-paths is only related to the grains size d (or compaction) and not to the tilting angle θ or the flow-rate Φ. This supposition could be in agreement with the results found for the diffusion coefficient D(θ) which is invariant under Φ and θ. x (cm) 0 5 -5 -10 10 - 15 15 x (cm) 0 5 -5 -10 10 - 15 15 x (cm) (a) At the beginning, the exploration is distributed along all the surface (due, in principle, to the 'percolation regime'). After some time, the exploration stabilizes in a well-defined zone (x ≃ [-9, 9] cm). (b) In order to obtain the local events, we calculate the difference between two consecutive images (t n+1 -t n ). The intensity indicates that the displacement goes from -1.5 cm to 1.5 cm (scale bar). (c) Binarizing this last result helps to observe the exact position of the local variations. The red color represents the emission of a bubble while the blue displays the static zones. We observe that the dynamics at the surface starts from different points (as a consequence of the initial 'percolation regime' injection) and in time concentrates around the center. The addition of all the values from the columns Σn z and rows Σn x , represent the total number of bubbles emitted at a given position and time respectively. The bubble emission distribution is centered around the inlet and the emission of bubble in time is constant after the stabilization of the air path [h g = 22 cm, h w = 2 cm, d =250-425 µm, Φ = 0.91 mL/s].

Discussion and conclusion

Results and applications

In this chapter, we have reported preliminary results concerning the formation, growth and finally the stationary state of the fluidized region. We have shown that the grains transport is particularly complex, observing a combination of motion (grains flow due to convective rolls and stick-slip). The associated horizontal speed of the interface evolves as a logarithm of time, and the data collapse when normalized by the square root of the height z. The dependence of D on θ is still difficult to interpret, and might be a direct consequence of the inner dynamics of the system, including memory effect and inertia. The role of the inertia has been only quantified for θ = 0 • and not for other angles. This may be important, specially when the gravity is smaller.

The dynamics inside the fluidized zone has shown two mechanisms that are not easy to explain. The first one, is a convective roll movements by the grains. The velocity, depends on the horizontal position, decreasing its magnitude when far to the center. The second one is the stick-slip zone (SSZ) observed at the boundary of the fluidized zone. In time the SSZ get thinner to finally desapears by the invasion of the fluidized zone. Interestingly, a similar description has been reported for eruptions in fluidized beds [START_REF] Gernon | Gas-fluidisation in an experimental tapered bed: Insights into processes in diverging volcanic conduits [Proceedings Paper[END_REF], with application to the formation of massive volcaniclastic kimberlite (MVK) (Brown, Field, et al., 2008). Figure 4.8 shows a schematic representation of MVK. As we can observe, the authors suggest, that after the fluidization of the conduit, the material (rock, and layered volcaniclastic kimberlite) focus toward the centre of the bed and then moves downward in a manner similar to that observed in our case.

Future improvements

For future work, we can list some improvements that can be done in the experiment:

• In order to identify the mechanism behind the fluidized zone, it could be interesting to insert into the bulk some particles that act like tracers. This could be useful to determine the associated velocity inside the fluidized zone and could help to improve the contrast between the fluidized zone and its boundary. In this framework, as the global interface of the fluidized zone is not easy to detect6 , we could improve the determination of its contour, and get more robust results.

• One important question is how the dynamics will change if the system is constrained in every directions. In particular, if the grain layer height is constrained. Indeed, we have shown that the fluidization zone starts at the non-restricted Brown, Field, et al., 2008). -The material ascends in the fluidized conduit.

Rock material moves downward along the margins in a stagnant unfluidized region (Gernon, Gilbertson, et al., 2008). After the eruption, avalanches bring the material back to the conduit center. Avalanches occur along the crater walls.

boundary (free surface of grains). If all the system presents the same fixed boundary conditions, we can wonder if the fluidized zone will be observed, and, if so, how it will depend on the initial compacity.

Increasing the water height h w

In all the experiments presented, the water height h w is a constant parameter that we have not changed. Indeed, we have limited a small water height h g ≃ 2cm, to avoid the transport of the grains in the water. Nevertheless, observing the dynamics of the free surface, we have noticed that the expulsion of air (from the bubble emission) ejects some grains that are then deposited back on the surface. The displaced grains slightly deform the surface around the vertical of the emission point. At short time scale, small craters appears, that are then destroyed after the expulsion of another bubble next to the previous one.

An interesting question is what happens if the water height h w becomes large enough for the grains to be freely advected in the water. In particular, how the surface shape will change if the advected grains are deposited not only back over the fluidized zone, but also on the outer sides. This question will be the central point of the next chapter.

In this chapter, we report the formation of a crater at the free surface of grains. when we increase the water height h w above the granular bed. We account for the dynamics of the crater formation in regard to the gas flow regimes. First, we describe the different steps that characterize the crater formation (Sec. 5.2.1). We then describe qualitatively the gas flow regimes (Sec. 5.2.2). Then, we describe the crater formation (Sec. 5.3) and make use of these observations to estimate the profile of the grain deposition. We also discuss the effects of the finite depth of the granular bed and of the finite water height (Sec. 5.3.3). We finally quantitatively study the influence of the gas flow-rate (Sec. 5.4.1) and of the grain size (Sec. 5.4.2).

Experimental Setup

In order to account for the dynamics of the free surface, we use the same experimental setup described in Sec. 2.1.2 [Fig. 2.3]. Additionally, we use a second configuration in order to observe the motion of the grains in the water. We light up the sample from the top with a linear light source (Polytec, DCR3) and take pictures with a high-resolution digital-camera (Nikon, D200). Initially, a thin layer of colored grains can be deposited at the free surface of the initially flat and horizontal bed.

Qualitative observations.

Before reporting quantitative measurements of the crater dynamics, let us mention some qualitative features of the crater formation.

Crater formation

The formation of the crater can be described as follows [Fig. 5.1]. First, air is injected at constant flow-rate at the bottom of the cell. After crossing all the immersed granular bed it reaches the grain-water free surface. Grains are lifted up by the ascending bubbles and the liquid flow which mainly consists of two large convective rolls produced by the bubbles. When the bubbles get to the water-air interface, the resulting flow (due to the convective rolls) pushes the grains away from the center. As the grains are subjected to gravity, they deposit back at a finite distance from the center onto the grain-water free surface. At this step, there are two possibilities for the grains. If they deposit on the inner side of the pile and if the local slope is larger than the angle of avalanche, the grains will flow downwards along the slope and will be lifted again by another bubble. On the contrary if the grains deposit on the outer side of the pile, they will contribute to the crater growth. The repetition of this process in time leads to the formation of two piles forming a crater.

Gas flow regimes and bubble emission mechanism

From the signal pressure and depending on the flow-rate, we obtain three qualitative regimes when the air crosses the granular bed, that we describe as follows:

Small flow-rate Φ -Bubbling regime • This regime is characterized by a regular emission of successive bubbles, independent from one another1 . The pressure signal exhibits successive rises and drops (sawtooth signal), the latter being associated with bubble emissions from the inlet at the base of the granular layer.

In the bubbling regime, a gas bubble, while growing underneath the free surface, pushes up a thin layer of grains which forms the bubble head. Once it has crossed In the bubbling regime, a large bubble forms underneath the free surface of the granular bed and then bursts. The explosion pushes out of the bed a large number of grains, which are subsequently advected upwards in the bubble wake. (c) In the open-channel regime, the walls of the channel are stable and a continuous gas flow escapes the granular bed. In this regime the grains are only torn out from the granular bed by the water flow and subsequently advected in the wake of the small bubbles that form in water [d = 400 µm, h g = 8 cm, h w = 16 cm, and Φ = 2.5 mL/s].

the interface, the bubble, while it rises up in the water, advects the grains in its wake [Fig. 5.2(b)]. As we can observe from the image, this regime is very efficient in the transport of grains.

Large flow-rate Φ -Open-channel regime

• The system sustains a continuous air-flow through a channel crossing the whole granular layer. The overpressure ∆P associated with this continuous air emission is almost constant.

In the open-channel regime, the air is released continuously through the channel whose walls remain at rest. The grains are advected upwards, from the free surface, by the water flow behind the small ascending gas bubbles [Fig. 5.2(c)].

In this case the transport of grains is less efficient.

Intermediate flow-rate Φ -Intermittence

• We observe a spontaneous alternation between the two previous regimes. The channel forms after the emission of several bubbles and subsequently pinches off after a finite time, leading the system back to the bubbling regime. As a consequence, activity and rest periods are observed in the pressure signal [Fig. 5.2(a)]. This phenomenon is explained by the ability of the material to sustain a stable channel thanks to its peculiar rheology2 .

We notice that the deformation of the free surface of the granular bed does not seem to alter the gas flow process: as pointed out in section 2.3.2, the path followed by the air in the immersed granular layer, and the subsequent gas emission, is sensitive to local conditions only and not to what lies upon it.

From the bubble emissions mechanisms observation, one could expect the bubbling regime to lead to a faster growth of the crater, because the explosive bursting of the bubble apparently lifts a larger quantity of grains. In order to investigate this point, we shall report the quantitative study of the crater growth in a large range of Φ (Sec. 5.4.1).

Geometrical description of the crater

The ejection of the grains from the free surface and the subsequent deposition of the granular ejecta lead to the rapid formation of a crater, which is formed by two granular piles symmetrically positioned on each flank. Due to the permanent grain transport, the crater grows: The piles height increases while they move symmetrically away from the center (Fig. 5.3).

Temporal evolution of the crater

We observe that the dynamics of the crater formation is very fast at the beginning and slows down through time. The results from the image analysis show that the distance L between the two summits is a logarithm of time which explains the fast growth at early stages [Fig. 5.4]. Figure 5.4 (a) display the space-time diagram of the free surface amplitude for one realization. The results indicate the symmetrie evolution of the two piles, that in some cases, can show slight differences in amplitude from one pile to another (the left pile is ∼ 0.5 cm higher than the right one). Figure 5.4 (b) shows the distance between the summits L(t) as a function of time t. From the inset we can observe the dependence of L ∼ ln(t).

Angle Evolution

Let us denote α and β the maximum angles that the piles make with the horizontal, respectively, inside and outside the crater [Fig. 5.3 and 5.5]. From direct observations, we notice that α and β can differ significantly, especially at the early stages of the crater formation [Fig. 5.5(b)]. We note that the inner flanks of the crater are continuously avalanching and exhibits a clear slope, and that α (about 28 • ) remains constant during the whole crater growth. In contrast, the outer slopes are not well-defined3 , especially at the early stages of the crater formation, and β evolves in time.

Inner flank: angle α

• At the center the grains are very locally pushed away from the surface by the upward liquid flow. The local slope almost immediately exceeds the angle of avalanche and grains flow downwards on both inner flanks to replace the granular material missing at the center. Thus, if a grain deposits onto the inner flank, it flows downwards toward the center, which again explains why α remains almost constant and only fluctuates between the angle of repose and the angle of avalanche [Fig. 5.5 (b)].

Outer flank evolution: angle β

• In contrast, at the early stages of the crater formation, the outer flanks only result from the deposition of the granular material. The local slope of the outer flank being everywhere smaller than the angle of avalanche, a grain, once deposited, does not move anymore. However, at a finite time t, β reaches the value of the angle of avalanche, which results in surface flows along the outer flank. As a consequence, β evolves in time [Fig. 5.5(b)] : Starting from a small value at the early stages of the crater formation, β increases to reach a constant value after a finite time. We note here that β then remains slightly smaller than α, which is probably explained by the fact that α significantly exceeds the angle of avalanche due to the continuous grain flow along the inner flanks.

Finite size effects

We investigate the condition for which the crater growth is limited by the finite depth of the granular bed, h g , and by the finite water height, h w . We also study the behavior of the convective rolls and the relation with the crater size.

Dependence on the grains height h g

The crater growth depends on the material that is ejected from the center of the cell. If the granular bed is small, then, the material is ejected until a given time after which there are no more grains at the center and the crater stops growing. Taking into account the angle of avalanche and the conservation of the granular material volume, we can estimate the maximum accessible value of L, to be:

L g max = (2 + √ 2)h g / tan α. (5.1)
We report experimental results obtained for h g = 8 cm so that L g max ∼ 50 cm, which is larger than the lateral width of the cell. The growth of the crater is therefore not limited by the finite depth of the granular bed in our experiments.

Dependence on the water height h w

In the same manner, the crater growth is possible due to the deposition of grains advected in the water layer. The water height h w therefore restrains the height of the piles which can obviously not exceed h w . We estimate, the maximum accessible value of L, to be:

L w max = 2(1 + √ 2)h w / tan α. (5.2)
However, we checked experimentally that L w max largely underestimates the finite waterheight effect: Reporting L at a given, large, time t = 30 min as a function of h w , we obtain experimentally that the growth of the crater is limited by the water height for h w < L/4 (Fig. 5.6). Interestingly, we observe, in addition, that L(30 min) does not depend on h w if h w > L/4 : The crater growth is not altered by the water height if this latter is large enough. We report experimental results obtained for h w = 16 cm so that L w max ∼ 64 cm is larger than the lateral width of the cell. The growth of the crater is therefore not limited by the finite water height in our experiments.

Evolution of the convective rolls

We observe that the vertical size of the large-scale convective rolls compares to h w , whereas their lateral size L r compares to L/2 [Fig. 5.7 (a)]. One thus would expect the typical size L r to depend on L and thus, on time. However, comparing the typical size 2 × L r of the convective rolls with the crater size L at different times shows no correlation [Fig. 5.7 (b)]. In fact, L r is constant in time. This led us to think that the typical vertical distance to be taken into account is probably the distance over which a grain follows the bubble in its wake, in other words, the maximum accessed altitude h w . We report the distance L between the piles as a function of time t, for different air flow rates Φ (Fig. 5.8). We observe that L increases logarithmically with t and we write:

L L 0 = ln Φt V 0 (5.3)
with V 0 being a typical volume. In the bubbling regime, when the gas emission at the free surface consists of the periodic emission of independent bubbles at the free surface, a simple argument can account for the proposed dependance on the total emitted gas volume, Φ • t, alone: The dimensions of the piles, length L or height, are expected to depend on the total number N g of grains displaced by the water flow. Thus, N g being proportional to the number of emitted bubbles, which itself, at a given time t, is proportional to Φ • t, we expect L to be a function of Φ • t alone, as observed experimentally.

We observe that the relation holds true even at large Φ when the system exhibits, almost always, the open-channel regime. Thus, contrary to the intuition, the small bubbles emitted by the open-channel are as efficient as the exploding bubble of the bubbling regime in lifting the grains. This observation is again in agreement with the fact that the grains are mainly lifted by the water flow and not expelled by the explosion of the bubbles. logarithmically with time, independently from the gas emission regimes. Indeed, at small Φ, the gas flow regime consists mainly in the emission of independent bubbles, whereas it consists mainly in a continuous emission of gas at large Φ. From the interpolation of the experimental data with L/L 0 = ln(Φt/V 0 ), we estimate L 0 = (3.5 ± 0.7) cm and V 0 = (0.4 ± 0.2) mL for d = 400 µm (black curves) whereas we estimate L 0 = (3.5 ± 1) cm and V 0 = (0.02 ± 0.01) mL for d = 100 µm (grey curve) [h g = 8 cm, h w = 16 cm].

Dependence on the grain size d

From the interpolation of the experimental data, one can estimate the typical length L 0 and the typical volume V 0 . One observes that the length L 0 does not depend significantly on the grain size (the slope in the semilog plot remains almost unchanged) and we estimate L 0 = (3.5 ± 1) cm for d = 400 µm and d = 100 µm. In contrast, we estimate V 0 ∼ 0.4 mL for d = 400 µm and a significantly smaller value V 0 ∼ 0.02 mL for d = 100 µm. This can be observed graphically in Fig. 5.9. The typical size L r of the convective rolls for both grain size is more or less the same, while the volume of grains ejected are drastically different. For the small grain size [Fig. 5.9 (up)], the amount of grains lifted is significantly larger than for the large grain size [Fig. 5.9 (down)]. If we interpret the volume V 0 as the gas volume necessary to lift or move a given quantity of grains, one can easily understand that V 0 decreases with d. Indeed, the advection process, whose associated force scales like d, competes with the buoyancy force which scales like d 3 . As a result, large grains are more difficult to lift than small ones and the gaseous volume necessary to move them away is larger4 .

Spatial distribution of the deposition

At this point, it is particularly interesting to make use of these observations to estimate the spatial distribution of the deposited grains. To do so, let us now consider the local height h(x, t) of the free surface at the distance x from the center at time t [Fig. 5.10(a)].

Close to the center, h(x, t) exhibits a complex behavior which results from both the deposition and the surface flows. Far from the center, h(x, t), which results from the deposition alone, evolves linearly in time. From this latter observation, we deduce that the advection is not significantly affected by the crater growth and that the distribution of the deposited grains, far away from the center, is almost constant in time. At an intermediate distance from the center, we observe a change in the temporal evolution of h(x, t) : At small times, the dynamics is only due to the deposition whereas later, when β reaches the avalanche angle, the local dynamics results from both the deposition and the surface flows. The volume, Q(x), of granular material deposited at the distance x from the center per unit distance and per unit time is proportional to the initial slope ∂h(x, t)/∂t far away from the center (Fig. 5.10).

Seeking for a simple mathematical description of the deposition flow, we guess that the grains, pushed away from the center by the liquid flow, have a negligible probability to deposit back at the center and we propose to interpolate

∂h(x, t) ∂t ∼ x L c exp (- x L c
), (5.4) which makes it possible to extract a characteristic length L c at which the grains are deposited away from the center. The data are successfully interpolated by a x Lc exp (-x Lc ) for large x, and we get at a rough estimate L c = (1.13 ± 0.04) cm [Φ = 2.93 mL/s, h g = 16 cm, h w = 8 cm].

Crater growth model

The dynamic of the crater growth is produced by the advection of the grains and their subsequent deposition. In order to build a growth model, we take into account two ingredients: the evolution of the deposited grains volume, and the distribution of the variation of the grains height.

From the previous section (5.4.3) we have found that the spatial evolution of the deposited flow can be written, for large distance x: .5) As the total grain volume is conserved, the volume of the ejected region must be the same as the volume of the deposited region [Fig. 5.11]. Assuming that the inner angle α is constant in time, which is already corroborated (Sec. 5.3.2), we find the three unknown variables:

Q(x) = a x L c exp - x L c . ( 5
                 x 1 = L 2 tan α x 2 = √ 2 2 L (1 + √ 2) x 3 = L 2 (1 + √ 2) (5.6)
As the volume ( 1) must equal to 2 × volume ( 2), a simple geometrical analysis gives

z x h w L/2 1 α x 2 x 3 x 1 z g 2 Figure 5
.11: Graphical representation of the dunes movement-As the grains volume is conserved in the system, the two triangles surface ( 1) and 2 × (2) must be the same. z g denotes the crater depth from the initial flat surface.

the volume of the displaced grains v g [Fig. 5.11]5 .

v g (α, L) = 1 8 tan α (1 + √ 2) 2 L 2 (5.7)
Noting that only the deposition of grains on the outer flank of the pile contributes to its growth, we can combine equations 5.5 and 5.7: (5.8) which leads to the equation governing the pile growth,

d v g dt = ∞ L/2 Q(x) dx,
1 4 tan α (1 + √ 2) 2 L dL dt = a ∞ L/2 x L c exp - x L c dx.
(5.9)

Resolving the integral, we get: .11) Finally, integrating and then deriving, we obtain that asymptotically, L increases logarithmically with time t according to: .12) with L 0 = 2L c and Φ/V 0 = a(1 + √ 2) 2 /tan α. Estimating roughly L 0 and L c (Sec. 5.3) for both grain size we find: d = 400 µm L 0 = 3.5 ± 0.7 cm L c = 1.13 ± 0.04 cm (5.13) (5.14) which are in fair agreement with the theoretical expectation L 0 = 2 L c .

1 4 tan α (1 + √ 2) 2 L dL dt = a L 2L c + 1 exp - L 2L c . (5.10) Provided that L ≫ L c (which is almost always satisfied experimentally), exp L L c dL = a 4 (1 + √ 2) 2 tan α dt. ( 5
L(t) ≃ L 0 ln Φt V 0 , ( 5
d = 100 µm L 0 = 3.5 ± 1 cm L c = 1.38 ± 0.04 cm

Discussion and conclusion

We reported experimental data on the dynamics of a crater growth at the free surface of an immersed granular bed crossed by an upward gas flow. We observed that, due to the peculiar transport of the grains in the wake of the rising bubbles, the typical size of the crater increases logarithmically with time. The dynamics is not altered by the gas flow regime and only depends on the overall gas flow rate and the typical size of the grains. However, even if we clearly established that the global dynamics of the crater is compatible with the spatial structure of the grain deposition around the gas-emission locus, we are still missing a complete modeling of the grain advection by the turbulent flow, which deserves to be further investigated.

Conclusion & Perspectives

In this thesis we have combined experimental, numerical and theoretical analyses to study the dynamics of air crossing an immersed granular bed. The results explain the dynamics from the initial injection to the emission of air at the free surface, including the bulk dynamics. We have found that this process can be accounted for by a diffusionlike model. The associated diffusion coefficient depends on an unique parameter χ which compares the width of the distribution of capillary pressure to the hydrostatic pressure over a grain size (Chap. 2).

We also have presented numerical simulations for the injection of a constant volume of air (Chap. 3). When the air paths are superimposed, the results show a well-defined morphology of the invasion zone. Interestingly, its geometry does not only depend on χ but also on the injected volume N. We have found an expression that relates the effective diffusion coefficient to the aspect ratio of the invaded region. Repeating the simulation for a cell with a bottom free boundary condition allows us to have the entire invasion geometry, and thus identify the range of parameters for which the process is diffusive or not.

In Chap. 4, we have presented preliminary experimental results for the dynamics at the interior of the bulk. Searching for a local description of the grain dynamics, we pointed out its complexity. Indeed, we observed phenomena such as convective rolls inside the fluidized zone and a stick-slip zone at its border. These last observations require further investigation.

Finally, we have characterized the morphology of the free surface of grains when the water height h w is high enough to allow the transport of grains in the water (Chap. 5). The formation of a crater is reported. We found the evolution of the angles as well as the distance L between the two piles in time. The results show a logarithmic dependence of L on the air-flow Φ and time t. We propose an analytical model that relates the volume of grains transported, and the deposition far from the center. The experimental results are well described by the model.

Perspectives

The experiments, numerical simulations as well as the theoretical models can be improved in several ways. The following points propose further investigation to complete

Conclusion & Perspectives

this work.

• Experiments: In order to characterize the formation of the fluidized zone, and its dependence on the boundary conditions, we propose to restrain the grains movement at the free surface. In particular, it would be interesting to observe under which conditions the fluidized zone appears. We can also change the granular material. One suggestion could be to introduce hydrophobic particles to observe the dynamics of air in this type of medium. Finally, a thorough investigation of the effect of viscosity of the host fluid would have to be performed.

• Numerics: The simulations presented in this work are based on very simple hypotheses. They do not account for the dynamics of the system, the inertia of the air flow, and even do not simulate the free surface of grains. Each one of these features could be improved in order to have a better numerical simulation of the system. In this sense, molecular dynamics simulations could be helpful to reproduce and understand the dynamics of the grains both inside the bulk (e.g. stick-slip, convection rolls) and at the surface.

• Theory: The analytical model for the diffusion coefficient could be improved, and a more detailed description could be possible to account for the dynamics of the stick-slip zone and the description of the convective dynamics inside the fluidized zone.

Additional observations

Finally it is interesting to report a phenomenon that has been observed while preparing the initial condition in the 2D experimental setup (Sec. 2.1.2). In order to prepare the initial condition, the experimental cell was completely filled with water, closed, and then flipped upside-down and back again to eliminate any bubble trapped inside the medium. During this last process we observed two interesting phenomena that we did not describe in the text but which would deserve further investigation. the grains gently settle down. Then, a compaction front ascends toward the surface. The front is characterized by an apparent constant velocity, and a small bump at the center. Note that the free surface altitude decreases as the compaction front advances.

• Rayleigh-Taylor Instability: After tilting the cell, the granular bed stands above the water layer. As grains start to fall down, we observed a Rayleigh-Taylor instability at the free surface of the grains, which evolves until the complete deposition of the grains (Fig. 5.12). From the figure we observe three steps.

At the beginning we have an initial flat surface, which is then destabilized by gravity. Granular fingering, mixing and sedimentation are observed. Note that the fingers density changes while sedimenting. This type of instability has been studied experimentally and numerically for a similar system where the fluid is either air [START_REF] Vinningland | Granular Rayleigh-Taylor instability: Experiments and simulations[END_REF][START_REF] Vinningland | Size invariance of the granular Rayleigh-Taylor instability[END_REF] or a water-glycerol mixture (Niebling et al., 2010a[START_REF] Niebling | Sedimentation instabilities: Impact of the fluid compressibility and viscosity[END_REF]. The major difference with our case is the larger dimension in our experiment (specially the gap space), which makes possible to observe the following additional process.

• Compaction Front: As the grains start to sediment, we observe a second front which corresponds to the compaction of the grain suspension. This front starts at the bottom of the cell and ascends toward the surface while the less compacted zone continues to sediment (Fig. 5.13). We can explain this observation by the fact that the system is not a perfect 2D cell, allowing the grains density in the gap to vary. As grains sediment, the friction between them (and the walls) could be responsible for the propagation of the compaction front.

The two above phenomena are only examples among the complex dynamics which can be further investigated in immersed granular material. The contribution proposed in this thesis only opens the way to the investigation of new physical phenomena not observed in dry matter.

stretched at the open end of a cylindrical cavity [START_REF] Vidal | Acoustic signal associated with the bursting of a soap film which initially closes an overpressurized cavityexperiment and theory[END_REF]. A constant air flow is injected inside and the film inflates until it breaks. The overpressure inside the cavity, the film shape and dynamics are reported for different fluids and inflation flow-rates.

A.2 Experimental setup and procedure

We draw a thin film at the open end of a plexiglas tube that is connected to a syringe [Fig. A.1]. The tube used is 10 mm diameter and 8 cm high. The film is a mixture of commercial dishwashing (La vaisselle, Casino) (50%), glycerine (25%) and water (25%). The initial thickness of the film cannot be controlled, but we measure the film thickness using an spectrometer (Ocean Optics, USB650 Spectrometer)1 . The collected reflectance spectra was estimated between 350 nm and 1000 nm using a scan rate of 25 milliseconds, with an integration time of t = 307.2 ms (0.3 s signal averaging time and 1 nm data interval). From the interference pattern obtained by the spectrometer we can find the film thickness in time. Indeed the film thickness is obtained by calculating the interference of multiple rays over a flat surface:

e = m 2D n n 2 -sin 2 (θ) (A.1)
where e is the film thickness, m the number of fringes, n the refractive index, θ the angle of incidence and D n the wavenumber region used (λ max -λ min ). Fig. A.2 (a), shows an example of an interference pattern at a given time. We detect five peaks in a range of 455 to 722 nm. The film has an initial thickness of ∼ 6 µm which decreases to ∼ 1 µm before bursting2 . As observed, the film remains for a long period of time at a thickness of about 1.2 µm. In the following, we will consider the film thickness e to be constant. Note that this hypothesis is rough, as the film thickness varies both in space [START_REF] Bergeron | Forces and structure in thin liquid soap films[END_REF] and time [Fig. A.2 (b)].

Protocol

Air is injected into the tube at a constant flow-rate controlled by a syringe pump ranging from 0.5 mL/min to 26 mL/min until the explosion of the film. The overpressure inside the tube is monitored by a differential pressure sensor (LPM5480 Druck) with a sensibility of ± 0.1 Pa. The film is observed with a webcam (Logitech, QuickCam Express, recorded at 20 fps) fixed on one side of the tube, which is illuminated by a lamp set behind the sample. This allows us to record the evolution of the bubble Air is injected in the system with a syringe and the bubble (blue line) grows. The arrows indicate the external forces applied to the system (pressure and surface tension). A zoom on the film shows the bubble interface. We consider a constant thickness e, R is the radius of curvature, r the cylinder radius, V the volume of the bubble, P 0 the atmospheric pressure and ∆P the overpressure inside the tube.

P 0 R r P 0 + ΔP θ V P 0 P i γ i γ e R
by detecting its edge. In order to compare the overpressure ∆P inside the bubble with the radius of curvature R (by detecting the edge bubble), we used a LED3 to synchronize the images and pressure acquisition: an electrical pulse, recorded together with the pressure signal, turns on and off the LED on the image sequences. We can then associate the radius of curvature of the bubble with the overpressure inside the tube at each time.

Usually for a high gas volume in the chamber the edge of the film (meniscus) in contact with the tube does not move significantly. To ensure that it remains fixed at this boundary, the tube was bevelled. Fig. A.3 shows an example of the film evolution. To obtain the radius of curvature of the film, we analyzed the images and detected the film evolution. We then fit the interface with a circle, assuming that this is the characteristic shape of the film evolution. The pressure signal is measured in time. We perform for each set of data three different series of measurement for the same flow-rate Φ, then repeat the process for all the flow-rates in order to get a good estimation of the error.

A.3 Dynamic pressure evolution

First, we compare the overpressure evolution inside the tube with the one predicted by the Laplace equation. Then, in order to account for the dynamics, we write a simple analytical model (Navier-Stokes) to quantify the dependence of the evolution of the A.4: Comparison between the theoretical and experimental overpressure -Comparing the experimental data (red circles) with the overpressure evolution predicted by Laplace equation (black line), we see that they collapse independently of the flow-rate as long as the volume of gas injected is small. Inset: for a large injected volume (typically 10 V 0 , where V 0 is the volume of the half sphere) the two curves split apart. In this case, the bubble is not spherical anymore (see text).

bubble pressure on the inertia, viscosity and capillary pressure. The main hypothesis in the following discussion is that we considerer the bubble thickness to be constant4 .

A.3.1 Comparison with Laplace pressure

We compare the pressure evolution inside the bubble with the overpressure predicted by the Laplace equation:

∆p = γ 1 R 1 + 1 R 2 (A.2)
where γ is the surface tension of the interface (in our case liquid-air) and R 1 , R 2 are the principal radii of curvature. To simplify the problem we consider a spherical bubble, i.e. the two radii of curvature being the same (R 1 = R 2 = R). As we have two interfaces (air-liquid inside the bubble, liquid-air outside the bubble), the capillary pressure can be written ∆p = 4γ/R. We can thus estimate the quasi-static evolution of the capillary pressure for the thin film.

We compare the static evolution of the pressure with the data (Fig. A.4). Both curves collapse, except for large volumes where the experimental data separates from the static evolution [Fig. A.4 (inset)]. The pressure difference in those regions can be explained by the fact that the film is not spherical anymore. Indeed when the volume of gas injected is large (typically higher than 10 V 0 , where V 0 = 2 3 πr 3 is the volume of a half of a sphere), the bubble is no longer spherical and changes to a bulb shape.

A.3.2 Analytical model: dynamic effects

In order to quantify the effects produced by the dynamics of the bubble, we calculate and contrast the dependence on the inertia, capillary and viscous terms. To find the dependence of the overpressure on the velocity and viscosity we use the Navier-Stokes equation for an incompressible Newtonian fluid:

ρ d v dt = -∇P + η∇ 2 v ∇ • v = 0 (A.3)
where v is the velocity of the bubble, ρ the liquid density and η its viscosity. We are interested in finding the pressure term P . In order to do so, we calculate each term of the equation. By assuming the liquid incompressibility, the fluid velocity v(r, t) can be written:

v(r, t) = R 2 r 2 Ṙ , (A.4)
where R(t) and Ṙ(t) are the bubble radius and velocity, respectively. First, we calculate the viscous term:

∇ 2 v = 1 r 2 ∂ ∂r r 2 ∂v ∂r = 1 r 2 ∂ ∂r r 2 R 2 Ṙ -2 r 3 = -2R 2 Ṙ r 2 ∂ ∂r 1 r = 2 r 4 R 2 Ṙ (A.5)
The next term is the unsteady acceleration:

∂ v ∂t = 1 r 2 2R Ṙ2 + R 2 R . (A.6)
Combining the two above terms, we can write the pressure gradient:

∇P = ρ r 2 2R Ṙ2 + R 2 R - 2η r 4 R 2 Ṙ (A.7)
Considering the radial symmetry, the integration gives:

P (r) = - ρ r 2R Ṙ2 + R 2 R + 2 3 η r 3 R 2 Ṙ + cte (A.
[N/m] η 10 -4 [Pa • s] e 10 -6 [m] R 10 -2 [m] Ṙ 10 -2 [m/s] ρ 10 3 [kg/m 3 ]
At last, using the Laplace term at the interior and exterior of the bubble ∆P i,e in equation A.8, we have:

P (r i ) = P i - γ R -e 2 = -ρ R -e 2 2R Ṙ2 + R 2 R + 2 3 η (R -e 2 ) 3 R 2 Ṙ + cte P (r e ) = P e + γ R + e 2 = -ρ R + e 2 2R Ṙ2 + R 2 R + 2 3 η (R + e 2 ) 3 R 2 Ṙ + cte (A.9)
Finally to obtain the pressure evolution we must find ∆P = ∆P (r i -r e ) (with the approximation e ≪ R, which is experimentally corroborated), we get: .10) where the first term corresponds to the Laplace pressure, while the second (viscosity) and third (kinetic energy) terms represent the dynamics effects. We can now compare each term in order to find which one is predominant. Using some general values for classic fluids (Table A.1) we find that the Laplace term is of the order ∼ 10 Pa, the viscosity term ∼ 10 -6 Pa and the kinetic term ∼ 10 -10 Pa. In order to observe any difference in our measurements we should therefore change the fluid viscosity over at least six decades or the flow-rate over ten decades.

∆P = 2γ R - ρe R 2 (2R Ṙ2 + R 2 R) + 2γ R 2 e Ṙ (A
In the following, we quantify the dependence of the overpressure ∆P on the flowrate Φ over one decade only5 . Contrary to the above estimation, we observe a dependence of the overpressure variation on the injected flow-rate Φ.

Dependence on the flow-rate Φ

In order to compare the pressure signals for different flow-rates Φ, we plot the overpressure evolution ∆P as a function of the normalized injected volume V nor = Φ • t /V 0 [Fig. A.5 (a)]. Interestingly, we find a difference in the amplitude of the signal. Note that the initial linear slope is the same for all the data, as expected, but the maximal pressure, ∆P = 4γ/r (where r is the tube radius) varies for the different curves. A.5 (b) shows the maximum pressure value ∆P max for each flow-rate. The results were calculated by taking the average of three series, with the error bar representing the standard deviation. The data reveal that there is a linear correlation between the maximal pressure and the flow-rate which could be the consequence of a dynamic surface tension. Note that the first fourth data points were taken with the same sample but another day. This could be the at the origin of the discrepancy with the global linear trend. results differ from what we expected (a same value), perhaps because the origin of the center of the bubble is not the same for all curves. Finally we should expect that γ evolves like the maximum pressure, meaning that there is a dynamic surface tension.

Radius of curvature and surface tension

A.4 Conclusion and perspectives

This complementary study revealed interesting informations on the physics of dynamic film bursting. We find that Laplace equation remains valid only by introducing a dynamic surface tension which is found to vary linearly on the injected flow-rate.

Additional measurements

• We performed additional measurements with a PDMS (Sigma-Aldrich, 200 ) of 10.000 cSt in order to compare the pressure evolution of another film with a high viscosity. Unfortunately, to draw a 'thin' film for high viscosity fluids is not an easy task to do. Furthermore, the results obtained were not reproducible and the bubble burst at any time. Indeed, the thickness of the film is not constant. From Fig.

A.6 we observe a snapshot for a high viscosity bubble. Note that the bubble thickness is not constant over the film, being thicker at the bottom and thinner at the top.

Many things would be interesting to pursue, starting with the improvement of the automatization of the pressure signal and the film detection. One important question arises from the maximum pressure detected for the different flow-rates. As we used materials that are not tabulated (typically dishwashing and glycerin) we do not have any control on the surfactants, which can induce some unknown behavior. It would therefore be necessary to repeat the experiment with well-controlled soap mixtures.

The analytical model could also be improved by considering a time-dependent film thickness. Finally, it would be interesting to correlate the pressure measurements inside the tube with acoustic measurement outside, as it has been done in the case of static bursting [START_REF] Vidal | Acoustic signal associated with the bursting of a soap film which initially closes an overpressurized cavityexperiment and theory[END_REF], in order to see if there is any correlation between the film thickness e and the acoustic energy at bursting. We report the formation of a crater at the free surface of an immersed granular bed, locally crossed by an ascending gas flow. In two dimensions, the crater consists of two piles which develop around the location of the gas emission. We observe that the typical size of the crater increases logarithmically with time, independently of the gas emission dynamics. We describe the related granular flows and give an account of the influence of the experimental parameters, especially of the grain size and of the gas flow. DOI : 10.1103/PhysRevE.79.021301 PACS number͑s͒: 83.80.Fg, 47.57.Gc, 47.85.Dh

I. INTRODUCTION

Craters are part of the widespread phenomena observed in nature, going from impact meteorite craters to volcanic structures. Studies of crater morphologies have a wide range of applications, going from puzzling crater formation in drying paint ͓1͔ to molecular dynamics ͓2-4͔. Among the main applications to natural phenomena, aside from meteorite impact crater, are the formation and growth of volcanic edifices, by successive ejecta emplacement and/or erosion. The time evolution and dynamics play a crucial role here, as the competition between volcanic-jet mass-flux ͑degassing and ejecta͒ and crater-size evolution may control directly the eruptive regime ͓5͔. Attempts have been made to model the talus development of volcanic caldera through erosion ͓6͔, or to constrain the morphology and dynamics of pyroclastic constructs via granular-heap drainage laboratory experiments ͓7͔.

Crater morphology in dry granular material has been extensively studied, both experimentally and theoretically ͓8-13͔. Most of these studies investigate the final, steady, crater shape resulting from the collision of solid bodies with the material surface and scaling laws are derived ͓14͔. Note however that some authors also reported experimental study of the crater formation dynamics, including growth and collapse after impact, in the "single impacting body" configuration ͓15,16͔. In a recent work, Wu et al. have extended these studies to a particles stream impacting a dry granular bed ͓17͔.

In immersed granular material, one reports craters generated by an underwater vortex ring ͓18͔, involving fluidized ejecta dynamics, or underwater impact craters generated by landslide ͓19͔. Craters in immersed granular materials can result either from two-phase or three-phase flows. In particular, water or gas flowing through an immersed granular bed can induce localized instabilities and fluidization ͓20,21͔, which eventually leads to the formation of craters at the free surface.

In a previous experimental study, Gostiaux et al. ͓22͔ have investigated the dynamics of air flowing through an immersed granular layer. They reported that, depending on the flow rate, the system exhibits two qualitatively different regimes: At small flow rate, the bubbling regime during which bubbles escape the granular layer quite independently one from another; at large flow-rate, the open-channel re-gime which corresponds to the formation of a channel crossing the whole thickness of the granular bed through which air escapes almost continuously. At intermediate flow rate, a spontaneous alternation between these two regimes is observed. Interestingly, they noticed the appearance of a crater around the locus of air release but did not provide any extensive study of its formation dynamics. Here, we extend these seminal observations to a quantitative study of the resulting deformation of the free surface of the granular bed. To do so, we reproduce the previous experimental conditions in two dimensions: In a vertical Hele-Shaw cell, the crater is then formed by two sand piles which grow and move away from each other as time passes by. By monitoring the evolution of the free surface through time, we investigate the effect of the different gas flow regimes on the crater dynamics. The results point out a grain-transport mechanism that differs significantly from the one involved in dry-sand dune motion ͓23-26͔.

II. EXPERIMENTAL SETUP AND PROTOCOL

The principle of the experiment is to force air, injected locally at the base of an immersed granular bed, to flow vertically across the material. The experimental setup thus consists in a 2-mm-thin vertical-cell containing the grains and water ͑Fig. 1͒. The two vertical walls ͑glass plates͒ are 40 cm wide and 30 cm high.

The air injection is assured by a system similar to that used in ͓22͔: Air is injected, at a constant flow rate, in a chamber connected to the bottom of the cell. In practice, a small hole ͑diameter 1 mm͒ connects, at the center, the bottom of the vertical cell to the chamber, partially filled with water in order to tune the inner gaseous volume V. The air flow is controlled thanks to a reducing valve which maintains a constant overpressure at the entrance of a capillary tube: Provided that the pressure drop in the capillary tube is much larger than the fluctuations of the pressure inside the chamber, the air flow toward the chamber is not significantly altered by the dynamics of the granular layer and remains constant. The flow rate ⌽ ͑ranging from 1.5 mL/ s to 4.3 mL/ s͒ is tuned by changing either the capillary tube ͑rough͒ or the pressure difference imposed by the reducing valve ͑fine͒ and is subsequently measured ͑to within 0. Instruments, PCI-6251͒, monitors the pressure difference, ⌬P ͑to within 10 Pa͒, between the chamber and the outside atmosphere.

In order to account for the dynamics of the free surface, we image the system from the side. In a first configuration, a transparency flat viewer ͑Just NormLicht, Classic Line͒ located behind the cell is used to achieve a homogeneous lighting of the whole bed which is imaged with the help of a Webcam ͑Logitech, QuickCam Express͒ connected to a PC. A small software ͑Astra Image Webcam Video Grabber͒ makes it possible to take one image of the system every 10 seconds and, thus, to record the dynamics during several hours ͑typically 24 hours͒. In a second configuration, in order to observe the motion of the grains, we light up the sample from top with a linear light source ͑Polytec, DCR3͒ and take pictures with a high-resolution digital-camera ͑Nikon, D200͒.

Initially, ⌽ is set to zero. Grains and water are introduced in the cell. The grains consist of spherical glass-beads ͑USF Matrasur, sodosilicate glass͒ that are sieved in order to control their size ͑typical diameters, 100-125 or 400-500 m͒. Initially, after sedimentation, the grains sit at the bottom and we make use of a thin rod to level the immersed granular layer. We denote h g the height of grains above the bottom of the container and h w , that of the water free surface above the granular bed.

The initial condition consists thus in a horizontal and flat layer of grains. After the opening of the valve, one observes that air is creating several paths between the hole at the bottom and the free surface of the granular bed. During this transient regime, some air channels merge, or cannot reach the upper layer of grains. ͑In this case, air bubbles remain trapped within the granular layer.͒ After several minutes, air crosses the granular bed along the vertical and bubbles are emitted in the water from an almost fixed position. Then, through time, a crater forms, the two piles on both sides growing and moving apart one from the other ͑see Fig. 1͒.

III. RESULTS

We shall report the dynamics of the crater formation in regard to the gas flow regimes. First, we shall describe quali-tatively the gas flow regimes ͑Sec. III A 1͒ and the associated mechanisms of the gas emission at the free surface ͑Sec. III A 2͒. Then, we shall describe qualitatively the crater formation ͑Sec. III A 3͒ and make use of these first observations to estimate the profile of the grain deposition. We shall also discuss the effects of the finite depth of the granular bed and of the finite water height ͑Sec. III A 4͒. We finally quantitatively study the influence of the gas flux ͑Sec. III B 1͒ and of the grain size ͑Sec. III B 2͒.

A. Qualitative observations

Gas flow regimes

As already reported ͓22͔, depending on the flow rate ⌽, two main regimes of the air flowing through the granular layer are observed. ͑The existence of the two regimes and the transition between them are tightly related to the non-Newtonian rheology of the granular material as proven by an experimental study conducted in another non-Newtonian material ͓27͔.͒ On the one hand, the bubbling regime, which is typically observed at small ⌽, is characterized by a regular emission of successive bubbles, independent from one another. In this regime, the pressure signal exhibits successive rises and drops, the latter being associated with bubble emissions from the injector at the base of the granular layer. On the other hand, in the open-channel regime, which is typically observed at large ⌽, the system sustains a continuous air-flow through a channel crossing the whole granular layer. The overpressure ⌬P associated with this continuous air emission is almost constant. In an intermediate range of ⌽, one observes a spontaneous alternation between the two regimes: The channel forms after the emission of several bubbles and subsequently spontaneously pinches off after a finite time, leading the system back to the bubbling regime. As a consequence, activity and rest periods are observed in the pressure signal ͓Fig. 2͑a͔͒. The phenomenon is explained by the ability of the material to sustain a stable channel thanks to its peculiar rheology. At last, we point out that the deformation of the free surface of the granular bed does not seem to alter the gas flow process: It has been previously FIG. 1. Experimental setup: Air is injected at constant flow rate ⌽ into a chamber ͑volume V͒ connected to a thin vertical cell containing an immersed, initially flat, granular layer. Through time, symmetrical piles form, grow, and move in opposite direction on either side of the crater. The experiment consists in observing the deformation of the granular-bed free-surface and in monitoring the overpressure ⌬P inside the chamber, in order to account for the influence of the degassing process on the crater dynamics ͓the picture here illustrates the bubbling regime ͑see text͔͒.

shown that the gas emission is mainly governed by local events occurring close to the free surface ͓22͔.

Bubble emission mechanisms

A close look at the free surface of the granular bed points out two different air-release mechanisms associated with the bubbling and open-channel regimes, respectively. In the bubbling regime, a gas bubble, while growing underneath the free surface, pushes up a thin layer of grains which forms the bubble head. Once it has crossed the interface, the bubble, while it rises up in the water, advects the grains in its wake ͓Fig. 2͑b͔͒. In the open-channel regime, the air is released continuously through the channel whose walls remain at rest. The grains are advected upwards, from the free surface, by the water flow behind the small ascending gas bubbles ͓Fig. 2͑c͔͒. From these observations, one could wonder about the relative efficiency of the two regimes in forming the crater. In particular, one could expect the bubbling regime to lead to a faster growth because the explosive bursting of the bubble apparently lifts a larger quantity of grains. In Sec. III B, in order to answer the question, we shall report a quantitative study of the crater growth in large range of ⌽.

Geometrical description of the crater

The ejection of the grains from the free surface and the subsequent deposition of the granular ejecta lead to the rapid formation of a crater, which is formed by two granular piles symmetrically positioned on each flank. Due to the permanent grain transport, the crater grows: The piles height increases while they move symmetrically away from the center ͑Fig. 3͒.

Before reporting quantitative measurements of the crater dynamics, let us mention some qualitative features of the crater formation. First, let us denote ␣ and ␤ the maximum angles that the piles make with the horizontal, respectively, inside and outside the crater. From direct observation ͑Fig. 3͒, one can notice immediately that ␣ and ␤ can differ significantly, especially at the early stages of the crater formation ͓Fig. 4͑a͔͒. We note that the inner flanks of the crater are almost straight and that ␣ ͑about 28°͒ remains constant during the whole crater growth. In contrast, the outer slopes are not straight, especially at the early stages of the crater formation, and ␤ evolves in time.

This behavior can be qualitatively accounted for by considering the grain flows ͑Fig. 5͒. Once lifted up by liquid flow, the grains are pushed away from the center by the liquid flow which mainly consists of two large convective rolls and then, subjected to gravity, deposit back onto the free surface at finite distance from the center. If the local slope is smaller than the angle of avalanche the grains do not move anymore once deposited. To the contrary, if the local slope is larger than the angle of avalanche, the grains flow downwards along the slope. At the center the grains are very locally torn off from the surface by the upward liquid flow. The local slope almost immediately exceeds the angle of avalanche and grains flow downwards on both inner flanks to replace the granular material missing at the center. Thus, if a grain deposits onto the inner flank, it flows downwards toward the center, which again explains why ␣ remains almost constant and only fluctuates between the angle of repose and the angle of avalanche ͓Fig. 4͑a͔͒. In contrast, at the early stages of the crater formation, the outer flanks only result FIG. 2. Overview on the gas flow process: ͑a͒ At intermediate flow-rate ⌽, the pressure signal exhibits a spontaneous alternation between activity and rest periods which correspond, respectively, to the bubbling and open-channel regimes. The emission of the bubbles at the free surface accordingly exhibits two qualitatively different behaviors. ͑b͒ In the bubbling regime, a large bubble forms underneath the free surface of the granular bed and then bursts. The explosion pushes out of the bed a large number of grains, which are subsequently advected upwards in the bubble wake. ͑c͒ In the openchannel regime, the walls of the channel are stable and a continuous gas flow escapes the granular bed. In this regime the grains are only torn out from the granular bed by the water flow and subsequently advected in the wake of the small bubbles that form in water ͑d = 400 m, h g = 8 cm, h w = 16 cm, and ⌽ = 2.5 mL/ s͒. FIG. 3. Temporal evolution of the crater: The black lines are the result of the automatic detection of the free-surface profile. The angles ␣ and ␤, denote, respectively, the maximum angle the inner and outer flanks make with the horizontal. We define L, the distance between the two peaks, and x, the distance from the center ͑d = 100 m, h g = 8 cm, h w = 16 cm, and ⌽ = 2.2 mL/ s͒. from the deposition of the granular material. The local slope of the outer flank being everywhere smaller than the angle of avalanche, a grain, once deposited, does not move anymore. However, at a finite time t, ␤ reaches the value of the angle of avalanche, which results in surface flows along the outer flank. As a consequence, ␤ evolves in time ͓Fig. 4͑a͔͒: Starting from a small value at the early stages of the crater formation, ␤ increases to reach a constant value after a finite time. We note here that ␤ then remains slightly smaller than ␣, which is probably explained by the fact that ␣ significantly exceeds the angle of avalanche due to the continuous grain flow along the inner flanks.

At this point, it is particularly interesting to make use of these observations to estimate the spatial distribution of the deposited grains. To do so, let us now consider the local height h͑x , t͒ of the free surface at the distance x from the center at time t ͓Fig. 4͑b͔͒. Close to the center, h͑x , t͒ exhibits a complex behavior which results from both the deposition and the surface flows. Far from the center, h͑x , t͒, which results from the deposition alone, evolves linearly in time. From this latter observation, we deduce that the advection is not significantly affected by the crater growth and that the distribution of the deposited grains, far away from the center, is almost constant in time. At an intermediate distance from the center, we observe a change in the temporal evolution of h͑x , t͒: At small times, the dynamics is only due to the deposition whereas later, when ␤ reaches the avalanche angle, the local dynamics results from both the deposition and the surface flows. The volume, Q͑x͒, of granular material deposited at the distance x from the center per unit distance and per unit time is proportional to the initial slope ץh͑x , t͒ / ץt far away from the center ͓Fig. 4͑c͔͒. Seeking for a simple mathematical description of the deposition flow, we guess that the grains, pushed away from the center by the liquid flow, have a negligible probability to deposit back at the center and we propose to interpolate ץh͑x , t͒ / ץt by a x L c exp͑-x L c ͒, which makes it possible to extract a characteristic length L c at which the grains are deposited away from the center. We shall show in Sec. III B that the shape of the distribution Q͑x͒ accounts for the evolution of the distance L͑t͒ between the two summits as a function of the time t ͓Fig. 4͑d͔͒. Temporal evolution of the crater geometry: ͑a͒ Inner and outer angles ␣ and ␤ vs time t ͑⌽ = 1.67 mL/ s͒. The grey zone marks the early stages of the crater growth during which the crater size compares with the size of the bubbles. The angle ␣, which is then not measurable, remains constant during the whole experimental time. In contrast, the angle ␤ increases continuously to reach a constant value after a finite time. In order to show that the slope is not constant along the outer flanks, we report the temporal evolution of the slope at x = We light up the sample from top with the linear light source and take pictures with the high-resolution digital camera. Initially, a thin layer of colored grains is deposited at the free surface of the initially flat and horizontal bed. Then, a series of ascending rolls pushes the grains away from the vertical central axis. Subsequently, the grains gently deposit back onto the free surface of the bed. Along the piles flanks, provided that the local angle exceeds the angle of avalanche, we observe continuous avalanches. Inside the crater, the flowing granular-material partly replaces the grains advected upwards at the center. Along the outer flanks, the deposited granular material either flows or sits at the free surface. One can clearly observe, in the inset, that new material, deposited far away from the center forms a thin layer of material which remains at rest. The deposition flow rate, Q͑x͒, is estimated from the temporal evolution of the growth velocity ץh / ץt of the deposited layer in that region ͓Fig. 4͑c͔͒ ͑d = 400 m, h g = 8 cm, h w = 16 cm, and ⌽ = 3.5 mL/ s͒.

Finite size effects

We aim at reporting the growth of the crater in an infinite surround, not limited by the finite depth of the granular bed, h g , or by the finite water height, h w .

The growth of the crater stops when the thickness of the granular bed at the center vanishes. Taking into account the angle of avalanche and the conservation of the granularmaterial volume, we can estimate the maximum accessible value of L, L max g = ͑2+ ͱ 2͒h g / tan ␣. We shall report experi- mental results obtained for h g = 8 cm so that L max ϳ 50 cm is about the lateral width of the cell. The growth of the crater shall thus not be limited by the finite depth of the granular bed.

In the same manner, the height of the piles can obviously not exceed h w . We estimate, the maximum accessible value of L, L max w =2͑1+ ͱ 2͒h w / tan ␣. However, we checked experi- mentally that L max w largely underestimates the finite waterheight effect: Reporting L at a given, large, time t = 30 min as a function of h w , we obtain experimentally that the growth of the crater is limited by the water height for h w Ͻ L / 4 ͑Fig. 6͒. Interestingly, we observe, in addition, that L͑30 min͒ does not depend on h w if h w Ͼ L / 4: The crater growth is not altered by the water height if the latter is large enough. We shall report experimental results obtained for h w =16 cm so that L max w ϳ 64 cm is larger than the lateral width of the cell. The growth of the crater shall thus not be limited by the finite water height. We report the distance L between the piles as a function of time t, for different air flow rates ⌽ ͑Fig. 7͒. We observe that L increases logarithmically with t and we write L / L 0 =ln͑⌽t / V 0 ͒, V 0 being a typical volume.

In the bubbling regime, when the gas emission consists of the periodic emission of independent bubbles at the free surface, a simple argument can account for the proposed dependence on the total emitted gas volume, ⌽t, alone: The dimensions of the piles, length L or height, are expected to depend on the total number N g of grains displaced by the water flow. Thus, N g being proportional to the number of emitted bubbles, which itself, at a given time t, is proportional to ⌽t, we expect L to be a function of ⌽t alone, as observed experimentally. In addition, we observe that the relation holds true even at large ⌽ when the system exhibits, almost always, the open-channel regime. Thus, contrary to the intuition, the small bubbles emitted by the open channel are as efficient as the exploding bubble of the bubbling regime in lifting the grains: This observation is again in agreement with the fact that the grains are mainly lifted by the water flow and not expelled by the explosion of the bubbles ͓22͔.

Dependence on the grain size d

From the interpolation of the experimental data, one can estimate the typical length L 0 and the typical volume V 0 . One observes that the length L 0 does not depend significantly on the grain size ͑the slope in the semilog plot remains almost unchanged͒ and we estimate L 0 = ͑3.5Ϯ 1͒ cm for d = 400 m and d = 100 m. In contrast, we estimate V 0 ϳ 0.4 mL for d = 400 m and a significantly smaller value V 0 ϳ 0.02 mL for d = 100 m. The precise dependence of L 0 and V 0 on d is difficult to access experimentally as such measurements require precise knowledge of the origin of time ͑a delay alters significantly the slope and the offset in the semilog plot͒ which we are missing. Nevertheless, V 0 is observed to decrease drastically with the grain size whereas L 0 remains almost constant.

IV. DISCUSSION

The crater growth is due to the advection of the grains and to the subsequent deposition away from the location of the gas emission. From the observation of the deposit far away from the center ͑Sec. III A 3͒, we propose that the deposited flow can be written, for large distance x, Q͑x͒ = a x L c exp͑-x L c ͒. This relation accounts for the logarithmic growth of the crater. Indeed, assuming the angle ␣ on The distance L between the summits increases logarithmically with time, independently from the gas emission regimes. Indeed, at small ⌽, the gas flow regime consists mainly in the emission of independent bubbles, whereas it consists mainly in a continuous emission of gas at large ⌽. From the interpolation of the experimental data with L / L 0 =ln͑⌽t / V 0 ͒, we estimate L 0 = ͑3.5Ϯ 0.7͒ cm and V 0 = ͑0.4Ϯ 0.2͒ mL for d = 400 m ͑black curves͒ whereas we estimate L 0 = ͑3.5Ϯ 1͒ cm and V 0 = ͑0.02Ϯ 0.01͒ mL for d = 100 m ͑gray curve͒ ͑h g = 8 cm, h w =16 cm͒.

both sides of the piles, a simple geometrical analysis gives the volume, v = 1 8 tan ␣ ͑1+ ͱ 2͒ 2 L 2 of the grains on the outer side of one pile ͑x ജ L / 2͒. Noting that only the deposition of grains on the outer flank of the pile contributes to its growth, we can write dv dt = ͐ L/2 ϱ Q͑x͒dx, which leads to the equation governing the pile growth,

1 4 tan ␣ ͑1 + ͱ 2͒ 2 L dL dt = a ͵ L/2 ϱ x L c expͩ- x L c ͪdx. ͑1͒
Provided that L ӷ L c ͑which is almost always satisfied ex-perimentally͒, we obtain that, asymptotically, L increases logarithmically with time t according to L͑t͒ӍL 0 ln͑⌽t / V 0 ͒ with L 0 =2L c and ⌽ / V 0 = a͑1+ ͱ 2͒ 2 / tan ␣. The experimental rough estimates, obtained for d = 400 m, L c = ͑1.13Ϯ 0.04͒ cm ͑Sec. III A 3͒, and L 0 = ͑3.5Ϯ 0.7͒ cm ͑Sec. III B 1͒ are in fair agreement with this expectation. The same conclusion holds true for the smaller grains, d = 100 m, as we measured L c = ͑1.38Ϯ 0.04͒ cm for the same value of L 0 .

Finally, if we interpret the volume V 0 as the gas volume necessary to lift or move a given quantity of grains, one can easily understand that V 0 decreases with d. Indeed, the advection process, whose associated force scales like d, competes with the buoyancy force which scales like d 3 . As a result, large grains are more difficult to lift than small ones and the gaseous volume necessary to move them away is larger.

The advection process remains difficult to model in details. Indeed, the liquid flow in the wake of the bubbles is generally turbulent and it is barely possible to account for the advection of the solid particles in such a complex stream field. There also remains a very important open question: We observed that the growth of the crater is not altered by the water depth ͑if large enough͒. However, we observe that the vertical size of the large-scale convective rolls compares to h w , whereas their lateral size compares to L / 2. One thus would expect the typical distance L c to depend on h w and, even, on time. Our study does not display such an effect. The typical vertical distance to be taken into account is probably the distance over which a grain follows the bubble in its wake, in other words, the maximum accessed altitude.

In conclusion, we reported experimental data on the dynamics of a crater growth at the free surface of an immersed granular crossed by an upward gas flow. We observed that, due to the peculiar transport of the grains in the wake of the rising bubbles, the typical size of the crater increases logarithmically with time. The dynamics is demonstrated not to be altered by the gas flow regime and only to depend on the overall gas flow rate and the typical size of the grains. However, even if we clearly established that the global dynamics of the crater is compatible with the spatial structure of the grain deposition around the gas-emission locus, we are still missing a complete modeling of the grain advection by the turbulent flow, which deserves to be further investigated. ACKNOWLEDGMENT G.V. acknowledges a grant by CONICYT ͑Comisión Nacional de Investigación Científica y Tecnológica, Gobierno de Chile͒.

I. INTRODUCTION

Fluid venting at the seafloor, usually concentrated along ocean margins, is a widely recognized phenomenon of geophysical, biophysical, and economical importance. On the one hand, the gas emitted is mainly hydrocarbon, whose extraction is of obvious economic interest [1,2]. On the other hand, the ecosystem on the seafloor close to the vents benefits greatly from the minerals brought by the gas or fluid emission [3]. Understanding the flowing regimes through a sediment layer is also fundamental from a geophysics point of view. Indeed, flow regimes in saturated granular media can go from homogeneous seepage (fluid percolation) to piping and partial failure [4]. This last case plays a major role in soil liquefaction, mud volcanism and diapirism, and hydraulic fracture [1,5].

The onset of piping phenomena, and its subsequent dynamics, has been studied extensively in the case of a two-phase flow (fluid and grains), both in physics and geophysics [4,[6][7][8]. It has been found that the channels are not stationary, but migrate across the surface, forming a complex subsurface plumbing system exhibiting different morphologies, representative of the fluid-venting dynamics [9]. At longer times, the channel migration fluidizes the region surrounding the loci of the gas emission, eventually leading to massive instabilities of water-saturated grounds [4].

Three-phase systems (solid, liquid, and gas) present an even more complex dynamics that has not been fully characterized up to now. For instance, in a two-dimensional (2D) porous medium that is not altered by the gas flow, one observes a transition in the morphology of the invasion front [10]. In granular systems, in which the grains are likely to move under the action of the gas flow (unconstrained immersed granular layer), three different regimes have been reported: a continuous bubbling regime during which one observes an almost periodic emission of gas bubbles at the free surface, an open-channel regime associated with a continuous gas emission through a stable channel that crosses the whole system, and an intermittent regime that corresponds to a spontaneous alternation of the degassing regime between the bubbling and the open-channel regimes [11][12][13]. Such dynamics is reminiscent of complex fluids, such as gels or polymers [14,15], crossed by ascending air flows.

In the present paper, we investigate the bubbling pattern that results from the gentle, punctual, injection of a gas at the base of an immersed granular layer. In this regime, the gas flow does not necessarily induce a significant motion of the grains. First, we characterize the gas-emission loci at the free surface of a cylindrical bed. Second, we report direct observations, in a 2D experiment, of the paths created by the gas in the bulk of the material. Finally, we show that a diffusionlike process accounts for the experimental observations and we relate the associated diffusion coefficient to the parameters of the problem.

II. SETUP AND PROTOCOL

The experiment consists in gently injecting air at the base of a column of immersed grains and observing the loci of the resulting bubble emissions at the free surface. The setup consists of a Plexiglas cylindrical tank (24 cm diameter) whose bottom is pierced at the center by an inlet hole (2 mm diameter) through which air is injected (Fig. 1). A reducing valve and a capillary tube are used to ensure that air is injected at a constant, small flow rate, , in a chamber of volume V , connected to the inlet (provided that the pressure drop in the capillary tube is much larger than the fluctuations of the pressure inside the chamber, the flow rate is not significantly altered by the dynamics of the granular layer and remains constant). The flow rate (ranging from 1 to 3 mL/s) is tuned by changing either the capillary tube (rough) or the pressure difference imposed by the reducing valve (fine) and is subsequently measured (to within 0.1 mL/s) by means of a homemade flow meter.

The initial condition is obtained as follows: the cylinder is filled with glass beads (USF Matrasur, sodosilicated glass) previously sieved in order to control their size (diameter d = 150-250, 250-425, and 425-600 μm). The height of the granular bed h g ranges from 2 to 24 cm. Distilled water is added and the material is vigorously stirred with a paddle to eliminate the gas bubbles trapped inside. The free surface of the granular bed is leveled and the water height above is h w ∼ 1 cm, which prevents the rising bubbles from hiding the emission locus.

Then, air injection is initiated by setting the flow rate to a chosen value. Air passes through the system and we observe at the surface a series of bubble emissions. To determine the loci of the gas emission, which vary in time, the free surface is imaged by means of a webcam (Logitech, QuickCam S7500) connected to a PC. The light source consists of a transparency flat viewer (Just NormLicht, Classic Line) placed slightly Air is injected at a constant flow rate into a chamber of volume V connected to the bottom of an immersed granular column (h g is the height of the granular column and h w is the height of the water above). The loci of the gas emission are determined by imaging the free surface by means of a webcam. out-of-axis to avoid direct reflections from the liquid surface into the camera. In such conditions, the emission of a bubble is marked by the appearance of a dark dot in the image of the free surface (Fig. 2, inset). A free software (Astra Image Webcam Video Grabber) is used to take 20 images per second with a resolution of 640 × 480 pixels. An automatic detection (Matlab, MathWorks) makes it possible to determine the loci of the gas emission with an accuracy of about 2 mm.

In the following sections, we report results obtained for 10-min acquisitions for different heights of grains h g , flow rate , and grain size d.

III. EXPERIMENTAL RESULTS

A. Qualitative description of the phase diagram

In a first step, we describe the different regimes observed in the experiment when varying the grain diameter d. We used three different grain sizes: diameters in the range 150-250, 250-425, and 425-600 μm.

For small grains (150-250 μm), when the air is injected, a large gas pocket forms below the surface. A dome growth due to the subsequent surface deformation is observed. Typicallyat any given air flow -the gas pocket grows for 3-5 min until the trapped air is suddenly liberated. In this case, the bubble emission at the surface occurs over a region limited by the size of the underlying gas pocket. The lower limit in the grain size can be understood as follows: when the grains are small, the Laplace overpressure γ /d is large enough to lift the whole granular column above as it exceeds the typical hydrostatic pressure ρ g gh g , where ρ g stands for the density of the grains.

For intermediate grain size (250-425 μm), bubbles are emitted at the surface of the granular layer, without any large-scale surface deformation. Bubble emission consists of successive, single bubbles, apparently randomly distributed in space around the centered air injection point (Fig. 2).

For large grains (425-600 μm), the behavior of the system is drastically different. The air crosses the immersed grain layer easily, producing not a single but multiple simultaneous bubble emissions at its surface. This transient behavior lasts for several seconds, up to several minutes, the air flow fluidizing rapidly the granular bed in a localized region until stabilizing to a single bubble emission location at the vertical of the inlet. This regime can be understood as follows: because of the large grain size, the capillary overpressure is small, and the air flows easily between the grains. The flow rate is not negligible anymore, and one cannot consider that the system is quasistatic and neglect the charge loss in the channels formed by the gas within the grains. The air then invades the material almost in the whole available free space between the grains.

This work focuses on the venting dynamics of single, successive bubbles emitted at the free surface of an immersed granular layer. Therefore, in the following sections, we will focus on grain diameters 250-425 μm only, and analyze the bubble distribution when varying the grain height h g .

B. Gas emission distribution

During 10 min, one observes several hundreds of events (Fig. 2). From the image analysis, we obtain the distance r separating the emission loci from the cylinder axis and report the corresponding values as a function of time t [Fig. 3(a)]. In the considered parameter range, we observe that events are widely, rather homogeneously, distributed around the vertical of the gas inlet. In addition, reporting the autocorrelation function R of the distance r [Fig. 3(b)], we observe that there does not exist any long-term correlation between the bubble emission loci. Indeed, R drops down quickly over a typical number of about two or three events, and, in any case, vanishes after 10 bubble emissions.

From the experimental data, we determine the probability distribution p(r) and the associated cumulative probability F (r) ≡ It is then particularly interesting to consider the dependence of ζ on the depth h g of the granular bed. For a given granular material, the data are compatible with ζ ∝ h g , as we will check in the following sections. This scaling does not depend 011302-2 significantly on the flow rate (Fig. 4). In the next section, we report the results of the 2D experiment, which makes it possible to visualize the air path through the granular layer and, thus, helps to understand the empirical law ζ ∝ h g .

C. Additional 2D experiment

In order to directly observe the paths of air within the granular material, we designed a two-dimensional setup that consists of a vertical Hele-Shaw cell: the immersed granular matter is contained between two vertical walls (glass plates 40 cm wide, 30 cm high, gap 2 mm). In the 2D experiment, the air flow is imposed by means of a mass-flow controller (Bronkhorst, Mass-Stream Series D-5111). The flow rate ranges from 0.2 to 2.2 mL/s. A webcam (Logitech, QuickCam Express) is used to image the system from the side, whereas the light source (transparency flat viewer, Just NormLicht, Classic Line) is positioned behind the cell. Again, the initial state of the system consists of an immersed granular bed, free of bubbles, whose surface is leveled. Air is then gently injected from an inlet at the center of the lower edge.

2D 3D ζ , ζ * (cm) 4 3 2 1 0 5 4 3 2 1 0 h g 1/2 , z 1/2 (cm
Direct observation from the side shows that air invades the system in two qualitatively different ways, depending on the flow rate. At small flow rate, the air does not induce any significant displacement of the grains while creating a new path and forms a rather complicated network of channels (Fig. 5). The granular material can then be, at least in the bulk, considered as a simple porous medium. However, due to local rearrangements that result from the emission of bubbles at the free surface, the channel network is permanently evolving. At a long time scale, the process leads to a region of smaller compaction at the vertical of the inlet. At larger flow rate (typically > 1 mL/s for d = 250-425 μm), one observes that the gas displaces the grains significantly while crossing the system (fingerlike behavior [16]). In this case, one usually observes the formation of a single channel that collapses after emission of the bubble at the free surface. The process repeats periodically, involving the formation of a different channel associated with each bubble.

As in the previous 3D experiment, we might determine where the bubbles are emitted at the free surface. In 2D, we chose a slightly different analysis. We suppose-and it is experimentally corroborated-that the path followed by the air in the immersed granular layer is sensitive to local conditions only (e.g., grain polydispersity and packing), and not to what lies upon it (grain height). Analyzing the air passing through a virtual horizontal line inside the 2D granular bed, located at height z, is therefore equivalent to considering the bubble emission distribution at the surface of an immersed granular layer of total height z. At long time, all the air paths create a fluidized zone (Fig. 6). Therefore, we chose to consider, in the 2D experiment, a deep immersed granular bed, and analyze the shape of the fluidized zone after injecting gas for several minutes or hours, until a steady state is reached. We can thus 011302-3 access, in a single experiment, to the characteristic size ζ * (z) of the region explored by the air as a function of the height z (Fig. 6). We observe that ζ * ∝ √ z, again independent of , which is compatible with the emission of bubbles at the free surface in a region having a typical size scaling like h g .

In the next section, we propose a theoretical description that accounts for the parabolic shape of the region explored by the bubble paths and, thus, for the dependence of the typical size ζ on the layer thickness h g .

IV. THEORETICAL APPROACH

In the present section, we report a simple 2D analytical model aimed at identifying the underlying mechanisms that govern the typical size ζ and show that the distribution of the bubble emission at the free surface can, as suggested by the experimental results, be accounted for by a diffusionlike process.

Qualitatively, let us consider first the propagation of air creating a path in the immersed granular material. A gas channel between the grains connects the path's upper end to the inlet at the bottom. Due to the injection at the bottom, the upper end moves when the inner overpressure overcomes the capillary overpressure associated with the narrow space between the surrounding grains. Depending on the local arrangement of the grains, the upper end can move upward or, more or less, sideward. Such a process results in a meandering of the path toward the free surface.

To describe the process more precisely, we consider, in a crude approach, that the air is creating a path in a 2D square array of grains [Figs. 7(a) and 7(b)]. To take into account the local rearrangement of the grains, we assume that the capillary overpressure associated with the possible directions of air path toward x < 0, δP -; x > 0, δP + ; and z > 0, δP z [Fig. 7(a)] are not identical but given according to a random distribution around a typical value δP c . In such a framework, potential correlations between successive bubble emissions are 1) and ( 4)].

neglected. Indeed, based on the experimental results reported in Fig. 3, which show that r is not correlated at long term, we can assume, in a first approximation, that each air channel forms in a medium that has not been altered by the previous one.

The system is subjected to gravity. Due to the hydrostatic equilibrium, there exists a typical pressure difference δP g between the surrounding liquid on top and on the sides. We can estimate δP g ρgd/2, where g stands for the acceleration due to gravity, ρ is the density of water, and d the typical grain size.

Thus, even if the capillary overpressures δP -, δP + , and δP z are equal on average, due to the hydrostatic overpressure δP g , the upper end of the path is more likely to propagate upward. The air will propagate upward if (δP z -δP g ) is smaller than both δP -and δP + . Let us denote by p z the associated probability. With such a definition, one can estimate that the upper end of the air path moves randomly N times to the left or to the right before moving once upward, with N ∼ 1/p z -1. As a consequence, the lateral mean-square displacement, associated with a displacement toward the free surface of the typical size d of a grain, can be estimated to be δx 2 ∼ Nd 2 . Thus, on average, the lateral distance ζ * between the path and the vertical of the inlet at the altitude z is given by ζ * ∼ √ D z, where To do so, let us consider that the capillary overpressures δP -, δP + , and δP z are distributed according to the distribution Q(δP ). For a given value of δP z , the probability for δP - and δP + to be larger than (δP z -δP g ) and, thus, for the air to propagate upward is [

D ∼ d 1 -p z p z . ( 1
∞ δP z -δP g Q(δP ) d(δP )] 2 .
Taking into account that δP z is also distributed according to the same distribution Q, we estimate

p z = ∞ -∞ Q(δP z ) ∞ δP z -δP g Q(δP ) d(δP ) 2 d(δP z ). (2)
In order to better interpret the general expression ( 2), one can consider the case of a Gaussian distribution:

Q(δP ) = 1 √ πσ p exp - δP -δP c σ p 2 . ( 3
)
The width σ p then accounts for the local disorder in the granular packing. In this case, we get

p z (χ ) = 1 π 3 2 ∞ -∞ e -v 2 ∞ v-χ -1 e -u 2 du 2 dv, (4)
where we defined χ ≡ σ p /δP g . Thus, p z and, as a consequence, D are functions of the unique parameter χ that compares the width of the distribution of the capillary pressure to the hydrostatic overpressure over the grain size.

The dependence of the coefficient D on dcan be easily evaluated numerically [Fig. 7(c)]: we observe that D decreases drastically when the grain size (respectively χ ) is increased (respectively decreased).

V. DISCUSSION

The simplistic model presented earlier indicates that the propagation of the air within the immersed granular bed can be accounted for by a diffusionlike model. The associated coefficient D = d f (χ ) is expected to scale like the grain size d and to depend on a unique control parameter, χ ≡ 2σ p /(ρgd), through an increasing function f (χ ). In this framework, one predicts ζ * ∝ √ z, compatible with ζ ∝ h g , in agreement with the experimental observations.

In addition, the diffusion coefficient is predicted to drastically depend on the grain size, not only because of the prefactor but also because of the dependence of χ on d. Indeed, χ = 2σ p /(ρgd) and, furthermore, σ p , which accounts for the distribution of the capillary overpressure, also depends on d: taking into account that the pore size (the space between the grains) is of the order of the grain size, one can estimate that the typical capillary overpressure δP c is of the order of γ /d, where γ stands for the surface free energy of the air-water interface. Denoting by w the typical relative variation of the pore size, one can estimate further that σ p ∼ w γ /d and, thus, that χ ∼ 2wγ /(ρgd 2 ). Thus, the control parameter χ and, in accordance, the coefficient D depend drastically on the grain size. For a granular bed made of monodisperse grains, the density can range, in practice, from that of the random loose packing, φ = 0.54, to that of the random close packing, φ = 0.64, so that w is of the order of a few percent. For a granular bed made of polydisperse grains, w is of the order of the relative width of the size distribution, about 25% in our typical experimental case in which d ranges from 250 to 425 μm. In this case, χ ∼ 10 and the model gives D ∼ 1.5 d. Even if we cannot expect from such a simplistic model a quantitative agreement with the experiments, we can evaluate that, for a typical depth of h g ∼ 25 cm, the region of the free surface in which the bubbles are expected to exit must be of the order of dh g ∼ 1 cm, thus centimetric, in fairly good agreement with the experimental observations (Fig. 4). We remark that, for χ > 1, D/d remains of the order of the unity and ζ ∼ dh g .

Due to the existence of very different regimes (Sec. III A), we have not been able to contrast our theoretical description with the experimental results in a large range of grain size. However, we do think that it remains valid for small grain size provided that the depth of the granular bed is large enough, so large that the capillary overpressure is not enough to lift the granular material as a whole.

The model we developed to describe the diffusionlike process is based on the hypothesis that all events (bubbles) are independent. There is, however, some correlation among the closest neighboring bubbles [Fig. 3(b)]. It points out the fact that a bubble rising in the granular medium creates a "weaker" path, that the next bubble will preferentially follow. This phenomenon has already been observed in complex fluids [17][18][19][20][21]. A way to improve the model could be by incorporating a weakness factor into the probabilities associated with the points in the mesh that were affected by the previous air path.

For large depth h g (typically, h g > 6 cm for d = 250-425 μm), we interestingly observe that, from time to time, stronger correlations between the loci of successive bubble emissions appear: the bubble emission remains in a very localized region, forming clusters of bubbles. Such events are revealed by small fluctuations of r around a finite value (Fig. 8). In this case, the air propagating through the immersed granular layer is able to create and maintain an open channel up to an intermediate height, the release of a bubble at the free surface being not enough to induce a collapse of the path along its whole length. A similar phenomenon (open channel) has already been reported, but the authors focused on the dynamics of the pressure signal in the chamber and not on the locus of the bubble emissions [11,12]. When a path partially bubble number n For large depth, a strong correlation between successive bubbles might be observed: several successive bubbles are emitted in a narrow region (cluster) before the system turns back to the regular regime (h g = 14 cm, d = 250-425 μm, = 4.3 mL/s.) 011302-5 stabilizes in the bulk, the upper end of the stable path plays the role of an intermediate inlet whose depth controls the size of the narrower region in which the bubbles are emitted at the free surface. After the collapse of the whole path, the system turns back to a more regular regime until, fortuitously, a new path partially stabilizes. Thus, for large depth, one might observe at the free surface clusters of emission loci instead of one single region centered at the vertical of the inlet.

VI. CONCLUSION

We showed that the spatial distribution of the bubble emissions at the free surface of an immersed granular bed can be accounted for by a diffusionlike process. As a result, the typical size of the region in which the bubbles are emitted is directly related to the depth of the granular bed h g and scales like D h g . The associated coefficient D scales like the typical grain size but is a rather complicated function of a unique parameter that compares the distribution of the typical capillary overpressure γ /d associated with the pores to the typical variation ρgd of the hydrostatic pressure difference over the grain size. As a result, the coefficient D is a decreasing function of the grain size: large grains lead to a narrow region at the free surface, whereas small grains tend to spread the loci of the bubble emissions. Finally, we interestingly observe that the partial stabilization of the path within the granular bed can lead to the appearance, at the free surface, of several clusters of bubble-emission loci, which might be relevant for explaining field observations [2].

The simplistic description of the model still deserves to be investigated further. To determine experimentally the relation between D and χ , taking into account that it is difficult in practice to characterize the granular bed (polydispersity and arrangement of the grains), we plan to vary, for a given granular material, the effect of gravity by tilting the 2D cell. The result of the latter study will be the subject of a future publication. We report a study of the paths formed by a finite volume of air gently injected at the base of an immersed granular material. A two-dimensional model, based on experimental observations, shows that the typical height and width of the region explored by the branched path depends not only on the injected volume V , but also on a dimensionless parameter χ which accounts for the relative effects of the gravity and capillarity. For a given injected volume V , larger gravity effects lead to taller and narrower structures; for a given χ , the typical height and width of the structure scale like V 1/2 and V 1/4 , respectively, while the typical gaseous fraction in the corresponding region increases accordingly like V 1/4 . Such results can be of practical importance: For instance, gas can be trapped on purpose in an underground natural container below a granular slurry. Our results can help in predicting if the gas is likely to reach the free surface and escape the system if the container presents a defect (hole or fracture The invasion of a gas in porous media is encountered in a wide range of systems, from industrial processes (oil industry [1,2], methane hydrate dissociation [3], etc.) to geophysical phenomena [4][5][6][7]. A typical example of the application is air sparging [8]: air is injected into the subsurface below the lowest known depth of contamination and, due to buoyancy, air serves to remove or helps degrade the contaminants. Concerning natural processes, the study of gas emission at the sea floor (venting dynamics) has led to the study of the release of methane from pockmarks in the mid-Atlantic continental shelf [9]. These studies contributed to the understanding of the carbon-cycle perturbations. The results are important because these pertubations are likely to induce global climate changes [10].

Morphology of air invasion in an immersed granular layer

Because of their implication in a wide range of systems and because of the potentially huge economical benefits, the understanding of such invasion processes has attracted the scientific community. They resemble many growth processes such as the Eden cluster model [11], ballistic models [12,13], dendritic growth, or the diffusion-limited aggregation (DLA) [14]. These models, which have been intensively studied [15,16], especially in numerical simulations, consider a homogeneous media in the absence of external forces [17].

One can also mention various numerical studies of similar systems, from the destabilization produced by gravity in two-dimensional (2D) porous media [18][19][20] to the study of the fractal dimension in an etched network [21]. Most of the systems involve tip effects in which the local growth velocity is proportional to the local gradient of an external field (pressure, impurity concentration, etc.). We also point out that a very similar growth or propagation phenomenon is observed in the mining industry when fragments are extracted from an extraction point at the base of the ore bed (draw body [22]).

Concerning, in particular, the injection of a gas in granular materials, we can distinguish the biphasic case (grains and air) [23] from the triphasic case (grains, fluid, and air) [24,25], the latter being more complex and not fully understood. The dynamics of the interface between air and an immersed granular material has been experimentally studied in a Hele-Shaw cell [26], mainly in a regime in which the gas forms a finger, reminiscent of the Saffman-Taylor finger. In this regime the grains are displaced by the interface. The main physical mechanisms at play involve the surface tension, the viscous dissipation, and a pressure gradient. The morphology produced during the air injection can be influenced by the rate of gas delivery and the vertical distance from the source [27]. We lately showed that, in a regime in which the air creates a path between the grains without moving them, the region invaded by the gas can be described as the result of a diffusion process [24]. The first simplistic approach neglected the formation of side branches during the growth of the air path. Here we propose a numerical analysis, based on experimental observations, of the role played by the injected volume and by the properties of the granular bed in the morphology of the region invaded by the gas, taking into account the formation of side branches.

II. PROBLEM STATEMENT

We aim at characterizing the geometry of the region invaded by a gas locally injected at the base of an immersed granular bed. In Ref. [24], we reported experimental results obtained in both three-dimensional (3D) and 2D experimental setups. The 2D system proved to be very useful as it makes it possible to visualize the paths created by the gas within the granular matrix. Thus, we will limit the present study to experimental examples obtained in a thin Hele-Shaw cell and numerical results from a (2D) model. We were previously interested in the loci of the gas emission at the free surface resulting from a continuous injection of gas. Here, we consider a qualitatively different situation in which a finite amount of gas is injected so that the invaded region does not reach the free surface. We aim at characterizing its morphology.

A. Summary and limitations of the former results

Long-term experiments [24] showed that, for the (2D) cell, the invasion of the immersed granular bed by the gas can be interpreted in terms of a diffusive-like model, the system exhibiting at long times a parabolic fluidized region 061302-1 1539-3755/2011/83(6)/061302 (7) ©2011 American Physical Society whose width w depends on the vertical distance z from the lower edge as w = √ D z. We proposed that the diffusion coefficient D was an increasing function of a unique parameter of the system χ ≡ σ P /ρgd, which compares the width σ P of the distribution of the capillary overpressure associated with the passage between the grains with the typical hydrostatic pressure variation ρgd over the grain size d. The result can be understood as follows: For very small grains, the hydrostatic pressure difference is very small compared to the capillary overpressure which scales like 1/d (thus much smaller than the width of its distribution) and, locally, the air creating its path between the grains propagates as if inside an isotropic medium. The width of the invaded region is large, which corresponds to large D. On the contrary, if the effect of the gravity is much larger than the capillary overpressure, the growth of an air finger in the system is not limited by the capillarity and the gas crosses the system straight along the vertical (very large grains). The invaded region is narrow, which corresponds to small D.

The previous experimental results are worth extending in several aspects. First, we were not able to check experimentally the dependence of D on χ . Indeed, one can barely vary χ in a well-controlled manner by changing the grain size d because the polydispersity of the samples and the packing of the grains are difficult to control. Second, in the theoretical approach, we assumed that the air was creating paths without side branches, which is obviously not the case when the effects of the gravity are weak. Finally, from a practical point of view, the geometry of the region invaded by a constant volume of gas is, at least, as interesting as the loci of the gas emission at the free surface in the continuous regime. Thus, we aim here at extending the previous results to branched invasion paths in the case of the injection of a finite volume.

B. Preliminary experimental results

To directly observe the paths of air within the granular material, we designed a 2D setup that consists of a Hele-Shaw cell (Fig. 1): the granular matter, immersed in water, is contained between two vertical walls (glass plates 40-cm wide, 30-cm high, gap 2 mm). To control the gravity effects, the cell can be tilted by an angle α with respect to the vertical so as to produce an effective gravity g eff ≡ g cos α. The experimental setup thus makes it possible to change χ , the granular material remaining unchanged. The injection of air through an inlet located at the center of the lower edge is insured by a mass-flow controller (Bronkhorst,. The flow rate can be tuned in the range 0.2 to 2.2 mL/s. The granular material consists of glass beads (USF Matrasur, sodosilicated glass) previously sieved to control their size (diameter d =150-250, 250-425, and 425-600 μm). A webcam (Logitech, QuickCam S7500, 640 × 480 px 2 , 20 images/s) is used to image the system from the side, whereas the light source consists of a transparency flat viewer (Just NormLicht, Classic Line) positioned behind the cell. The initial state of the system consists of an immersed granular bed (typical depth 24 cm), free of bubbles, whose surface is leveled. A chosen volume of air V is then gently injected in the system. The effective gravity is changed by tilting the cell by an angle α (from 0 • to 90 • ). The formation of the pattern in the granular bed is recorded with a camera positioned in front of the experimental cell.

The preliminary experimental results (Fig. 2) show that the branches are more numerous and, thus, the pattern more compact when the effective gravity is reduced by tilting the experimental cell. Accordingly, the invaded region is broader and smaller in this case whereas an increase of g eff leads to a narrower and taller pattern with less branches. To characterize the average geometry of the invaded region we would have to repeat the experiment a large number of times, which is difficult. Indeed, once the pattern is formed, to reset the initial condition, one must open the cell to remove the air trapped in the system, which takes a long time. Moreover, we would not be sure to prepare the system in the exact same way and the state of the granular packing (its density, for instance) might be different from one run to another.

Thus, to overcome the difficulty, we performed the numerical analysis of the problem that is thoroughly described in Sec. III.

III. NUMERICAL ANALYSIS

To assess the dependence of the geometry of the invaded region on the problem parameters, especially χ and the injected volume V , we perform the simple numerical analysis whose ingredients and results are reported in the Secs. III A and III B, respectively.

A. Ingredients

In a simplified 2D approach of the system, we consider that the granular packing reduces to a 2D square network, each vertex corresponding to the void space between four grains. The size of the network (1200 × 1000) compares with the typical size of the experimental cell in units of grains (typically 1600 grains in width by 1200 grains in height) and is large enough to avoid boundary effects on the sides and at the top. The only boundary condition is that the gas cannot trespass on the bottom plane.

The propagation of an air finger in the material is limited by the capillary overpressure δP c to overcome going from one vertex to a neighboring one. Considering that the typical size of the pass between the two void spaces scales like the grain diameter d [25], we estimate that, on average, δP c ∼ γ /d, where γ denotes the air-water surface tension. However, due to the local heterogeneity of the system (polydispersity, wetting conditions, local arrangement of the grains, etc.), δP c differs from one pass to another. Thus, to account for the heterogeneity, we assume that the links between the vertices are associated with capillary overpressures distributed according to a Gaussian distribution of width σ P around the nominal value P c which, we remind, is of the order of γ /d.

The propagation of the air finger along the vertical is favored by the additional contribution of the hydrostatics. Indeed, in the experiments, the water that fills the space between the grains is subjected to gravity so that, for instance, the pressure difference in the water between one vertex and the first neighbor above is δP g = ρg eff d, where ρ denotes the density of water. In the numerical computations, the effects of the gravity are accounted for by considering that the threshold overpressure to overcome to go from one vertex to another is δP t = δP cρg eff z, where z (positive) denotes the vertical distance from the outlet plane (bottom edge). Technically, we associate with each of the links the threshold overpressure δP t by adding the corresponding contribution of the gravity to the previously chosen map δP c .

The air path is calculated as follows: First, the finger is grown from the virtual outlet, at the center of the bottom edge. Second, during the finger growth, from an already existing path, we consider the whole set of links connected to vertices occupied by the gas and determine the one corresponding to the smallest value of the threshold overpressure δP t . We make the air invade the corresponding vertex, by adding the latter to the air path. The procedure thus allows the generation of branches (Fig. 3). We repeat the process until a chosen number of vertices (N = 2000,4000, and 8000), corresponding to a chosen volume V ≡ Nv p (v p is the typical volume associated with the void space between the grains), are occupied. Finally, to get a relevant estimate of the average geometry of the invaded region, we repeat the whole process 2000 times for the same parameters (for a given χ and N ).

B. Results

In the present section, we discuss the geometrical properties of the region invaded by the gas. The discussion is based on averages of 2000 numerical paths (Fig. 4).

Geometrical characteristics of the invaded region

As expected, for χ = ∞ [Fig. 4(a)], the gas invades the granular packing in an almost isotropic manner, the paths filling half a disk above the outlet. The only anisotropy results from the boundary condition (no air flow) at the lower edge of the system. The observed pattern resulting from a pure diffusive process with a reflecting boundary the typical size of the invaded region scales like d √ N . When χ is decreased, and thus the pressure gradient increased, the invaded region elongates along the vertical. However, one can notice that, for a given injected volume, the width of the structure does not change significantly whereas its height significantly increases ].

To report a qualitative behavior of the system, we first determine the contours of the invaded region associated with the isodensity lines in Fig. 4: A contour is defined by the fraction f such that, over 2000 paths, the contour is locally crossed 2000f times [Fig. 5 A quantitative characterization of the structure geometry is provided by measurements of its width x and height z. From the superposition of 2000 paths [Fig. 5(a)], we consider the intensity profile along the x axis, estimated over the whole height of the system [Fig. 5(b)]. As expected, the profile is symmetric and we consider x to be the width at half the peak value. In the same way, we consider the intensity profile along the z axis, estimated over the whole width of the system [Fig. 5(c)]. We observe that the profile is almost flat and suddenly decreases above a given altitude. We define z as the altitude of the point corresponding to half the plateau value.

Dependence on V and χ

Reporting the contours of the invaded region for various values of the injected volume (various N ) and effective gravity (various χ ), one observes that, qualitatively, the geometry of the pattern highly depends on both control parameters. First, one observes that, for increasing injected volumes [Fig. 6(a)], the height z increases almost linearly with N . We observe that z increases faster than the average width x, which results in an increase of the aspect ratio [Fig. 5(d)]. Second, for a given injected volume [Fig. 6(b)], a decrease in χ results in an increase in the height z and in a decrease in the width

x, which result in a drastic increase in . Note that, when χ is decreased the shape of the invaded region goes from a half-disk to a parabola or cone. We expect that the air goes straight up through the system for χ 1. The quantitative changes in the shape of the invaded region can be assessed by reporting as a function of χ [Fig. 5(d)]. For the whole set of experimental data, we observe that 5,50,166,500,1666, and ∞. The system exhibits a diffusive behavior for small χ . The curvature of the invaded region near the outlet is obtained from the slope at z = 0 (N = 8000).

∝ 1/

√ χ in a wide range of χ and reaches a plateau value of the order of the unity 0.82 for large χ . For a compact structure (isodensity), one would expect the aspect ratio to tend to 1 in the absence of gravity, in the absence of a lower boundary. Taking into account that the density is not constant within the structure and that the air flow is limited by the lower edge, we are not thus surprised that the ratio tends to a value of the order of the unity, but not exactly to 1. Interestingly, we note a collapse of the measurements when reporting the as a function of χ/ √ N : the scaling ∝ 1/ √ χ holds true for χ 30 √ N whereas 0.82 for larger values.

Additional pieces of information about the invasion process can be obtained by considering in Fig. 7(a) the probability density along x for a given altitude z: From the superposition of 2000 paths, one obtains the number of passages at a distance x from the axis, for a given altitude z, by considering the corresponding density profiles [Fig. 7(b)]. From the profiles, one gets the typical width σ of the region crossed by the air paths as a function of z [Fig. 7(c)]. We observe that, for small enough values of χ and z, σ 2 is proportional to z, which accounts for the parabolic shape of the invaded region near the outlet. By extension, whatever the value of χ , the shape of the invasion pattern near the outlet can be accounted for by its radius of curvature, equivalent to an effective diffusion coefficient D eff ≡ ∂σ 2 /∂z| z=0 [Fig. 8(a)]. We obtain

IV. DISCUSSION AND CONCLUSION

Contrary to classical growth models [11][12][13][14], our model takes into account both a heterogeneous medium (capillary overpressure distribution) and the effect of an external field (gravity). Previous studies of fluid invasion in a 2D porous medium, including the gravity destabilizing effect, focused on the geometry of a single invasion pattern [18,19]. By contrast, our work aims at characterizing the morphology of the region potentially explored by the invading fluid. The numerical results show that, as expected, the shape of the invaded region depends drastically on χ (i.e., on the effective gravity). Less obvious, an increase of N (i.e., of the injected volume) does not simply lead to a dilation of the invaded region, but to a change in (i.e., in the geometry).

The present numerical model exhibits a much richer behavior than that exhibited by our previous analytical analysis [24]. Indeed, formerly, the side branching and the 061302-5 return of the air path to a previous lower position were not allowed. The path was forced to grow upward. The main difference in the results is that, formerly, the radius of curvature (or effective diffusion coefficient) D eff was predicted to be a complex function of χ whereas the prediction of the present numerical results is that D eff d χ for χ 30 √ N . In addition, we report that the vertical size of the structure is linear as a function of the injected volume V and we estimate the dependence of the associated velocity on the control parameter χ , ν z 0.8d/ √ χ.

It is particularly interesting to interpret the parameter χ ≡ σ P /ρgd, which compares the width of the distribution of the capillary overpressure within the pores to the variation of the hydrostatic pressure over the grain size. On the one hand, for a given injected volume N , the air inside the existing path being connected, one can consider the path as an isobar. On the other hand, the pressure in the liquid, outside the air path, increases from the path tip (its highest point) downward, which makes the formation of the side branches more and more difficult at depth, below the path tip. Considering the meaning of χ , one can estimate that side branches cannot form at a distance larger than d χ below the tip. As a consequence, the gas injection results either in the formation of side branches in a region of typical height d χ (in a local increase of the gas fraction or in the widening of the structure) or in the growth of the tip upward. The result of such complex dynamics is a complex path exhibiting more or less side branches depending on the value of χ . Interestingly, on average, the gas occupies an elongated region whose radius of curvature near the outlet is d χ, as shown by the numerical results.

From the dependence of the height z and typical width √ D eff z, we can estimate the typical gaseous fraction F inside the invaded region. Estimating the corresponding surface area √ D eff h 3/2 , one gets from simple algebra F ∼ χ 1/4 / √ N , thus dependent on χ and on the injected volume (on N). The fraction F slightly increases when χ increases (i.e., when the gravity effects are reduced and the side branching enhanced). In addition, denoting w the typical relative variation of the pore size as proposed in Ref. [24], one can estimate further that σ p ∼ wγ /d and, thus, that χ ∼ wγ /(ρg eff d 2 ). Thus, considering that the pore volume v p ∝ d 2 and taking into account the result obtained for ν z , we are taught that the maximum height reached by the gas within the granular does not depend on the grain size and scales like V /(l c √ w) where l c ≡ √ γ /ρg eff . Thus, provided that the proposed estimate of χ is correct, for a given volume V , the maximum height is controlled by the capillary length l c and the relative width w, which account for the heterogeneity of the capillary overpressure.

In conclusion, we reported results of a numerical study which makes it possible to predict from the knowledge of one single control parameter χ the typical height, width, and gaseous fraction of the region invaded by a given volume of gas liberated at the base of an immersed granular bed. Such results could be of practical importance: For instance, gas can be trapped on purpose in an underground natural container below a granular slurry. Our results can help in predicting if the gas is likely to reach the free surface and escape the system if the container presents a defect (hole or fracture). The present study will be extended, from the theoretical point of view, to a slightly different geometrical situation, especially to the 3D case and, from the experimental point of view, to the case of a horizontal liquid flow.

A. 2

 2 Interference fringes and bubble thickness evolution . . . . . . . . . . . . A.3 Film bulging when injecting a constant air flow-rate . . . . . . . . . . . A.4 Theoretic and experimental comparison . . . . . . . . . . . . . . . . . . A.5 Radius of curvature and surface tension . . . . . . . . . . . . . . . . . . A.6 High viscosity bubble . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Tables 1.1 Grain size distribution analysis for the different batchs . . . . . . . . . 2.1 Parameters in the experiments and their associated dimension . . . . . A.1 Typical values for a soap film at room temperature. . . . . . . . . . . .

1 2Figure 1 . 1 :

 111 Figure 1.1: Overview model of fluid migration in the subsurface, from Naudts et al. (2008) -Carbonate formation with seep relocation and the acoustic sea-floor backscatter expression. Gas-bubble release from black venting holes. Note that the emission can be dispersed in a wide area or localized in a specific region.

Figure 1 . 3 :

 13 Figure 1.3: Examples of porous networks -(a) Rigid porous network. The matrix

Figure 1 . 5 :

 15 Figure 1.5: Models of digitations -A Porous medium can be described by a network where each vertex represents a pore. Depending on the property of the invading and host fluid [e.g. wetting (represented in blue) or non-wetting (represented in white)] we can have two types of invasion. When the non-wetting fluid displaced the wetting, the process is called drainage. The inverse process is called imbibition.

Introduction

Figure 1 . 6 :

 16 Figure 1.6: Experiments and numerical simulations in porous media (from Lenormand & Zarcone, 1989) -(a) Experiment: displacement of a wetting fluid (e.g. oil, black) by a nonwetting fluid (e.g. air, white). (b) Numerical simulation of a similar system (shown at the percolation threshold). The patterns resemble the experimental observation.

Figure 1 . 7 :

 17 Figure 1.7: Flow regimes in water saturated media (from Mörz et al., 2007) -(A) Homogenous Darcy type flow with water percolation bound to the fixed pore space. (B) Formation of micro-channels and conduits. Small-scale heterogeneities develop. The sample is still stable and macroscopic flow rates obey Darcy law. (C) Hydraulic failure of the sample. Pipe formation leads to a partial liquefaction.

Figure 1 . 8 :

 18 Figure 1.8: Homogeneous injection of air in dry granular material (from Nermoen, Raufaste, et al., 2010) -Depending on the injection flow-rate, different regimes are observed: (A) Darcy flow, (B) Nucleation of pipes and (C) full separation of the two media.

Figure 1 . 9 :

 19 Figure 1.9: Fluid emission location (from Zoueshtiagh & Merlen, 2007) -Photographs of water emission on the beaches of Georgiopoli in the Creta island (see text).

FracturingFigure 1 . 10 :

 110 Figure 1.10: Phase space [from Sandnes et al. (2011)]-Tentative phase diagram of morphologies in the ϕ -1 -Φ plane extended to the ϕ = 1 limit corresponding to a closepacked porous medium. Various displacement morphologies in the 'Frictional fluid' give way to fracturing ('Deformable media'), followed by a transition to capillary/viscous fingering ('Porous media').The phase boundaries are 'guides to the eye'.

Figure 1 . 11 :

 111 Figure 1.11: Pattern formation in a confined granular suspension (from Chevalier et al., 2009) -Increasing the grain fraction the resulting pattern formed by the air can transit from: (a) the classical Saffman-Taylor instability, to (b-c) pattern formation in suspension and finally (d) viscous fingering in a classical porous medium (see text).

Figure 1 . 12 :

 112 Figure 1.12: Particle size distribution of the glass beads used in our experiments -Each row show an image sample of the glass beads and the respective size distribution for: (up) d =200-250 µm, (center ) d =250-425 µm and (down) d =400-500 µm. The three grain size distributions are well-fitted by a normal distribution. We defined the average grain size as d = µ ± σ.

Figure 2 . 1 :

 21 Figure 2.1: 3D experimental setup -(a) Air is injected at a constant flow rate Φ at the bottom of the immersed granular column (h g is the height of grains and h w the height of water above the granular bed). The emission of bubbles at the free surface is recorded by a webcam positioned over the experimental cell. (b) Image of a bubble emission at the free surface [d = 250 -425 µm, Φ = 2 mL/s, h g = 12 cm].

Figure 2 . 2 :

 22 Figure 2.2: Typical overpressure signal -Pressure evolution for the bubbling regime. The emission of bubbles is characterized by a sawtooth wave [d =250-425 µm, Φ = 2.5 mL/s].

Figure 2 . 3 :

 23 Figure 2.3: 2D experimental setup -The 2D experimental setup makes it possible to observe the dynamics of air crossing the media. The air is injected at constant air flow (see text). In the experiment, the cell can be tilted from θ = 0 • to θ = 90 • , thus changing the effective gravity in the system. (b) Example of air crossing the granular bed [d = 250-425µm, Φ = 2.17 mL/s, θ = 0 • ].

Figure 2 . 4 :

 24 Figure 2.4: Gas emission detection -The white circles indicate the successive bubble emission at the surface of the immersed granular bed. The black circle indicates the typical size of the region explored by the air at the surface, centered at the vertical of the injection point at the surface (black cross). Inset: Example of a bubble emission image [d =250-425 µm, Φ = 3 mL/s, h g = 20 cm].

Figure 2 . 5 :

 25 Figure 2.5: Bubbles emission location as a function of time -(a) Distance of the emitted bubble from the center, as a function of time. (b) Cumulative probability F (1-F is represented, see text) and gaussian fit of the distribution of gas emission location (distance r from the center). The full width at half height of the Gaussian distribution is 2ζ, with ζ ∼ 2.02 cm [represented in (a) by the dashed lines]. Inset: Log-lin plot of the data (circles) and Gaussian fit (full line). (c) Autocorrelation R of the signal as a function of the neighboring bubble number n [d =250-425 µm, Φ = 3 mL/s, h g = 20 cm].

Figure 2 . 7 :

 27 Figure 2.7: Different invasion regimes at short time scale -(a) Percolation regime, usually observed at low air flow-rate. The air penetrates the medium as it would invade a rigid network. Note that we have simultaneous emission points at the free surface. (b) Fracture-like regime, usually observed for high air flow-rate. This regime is characterized by the displacement of grains causing the fracture of the medium. In this regime, only one single branch reaches the surface [Φ = 2.17 mL/s, h g = 22 cm, d = 250 -425µm].

Figure 2 . 8 :

 28 Figure 2.8: Dynamics of the crossing the 2D granular layer -(a) At short times scale, the air creates several paths within the granular bed that reach the granular free surface (here, percolation regime). (b) At long time scale, we observe the generation of a fluidized region, characterized by a lighter contrast produced by lower grain compaction. Its shape is determined by image analysis, and subsequently interpolated by a parabolic function (white dashed line). We denote by 2ζ * (z) the characteristic width, in the bulk, of the region explored by the air in the stationary regime at height z [h g = 22 cm, h w = 2 cm, d =250-425 µm, Φ = 0.43 mL/s].

Figure 2 . 9 :

 29 Figure 2.9: Model description -(a) Sketch of the arrangement of grains in the model. The air injected at the bottom propagates through the easiest path between the grains (dashed line). (b) Definition of the capillary overpressures δPfor each directions of the air path, x < 0 (δP -), x > 0 (δP + ) and z > 0 (δP z ). The probability for the air to propagate downwards is prohibited in the model.

Figure 2 . 11 :

 211 Figure 2.11: Normalized diffusion coefficient D/d as a function of χ -(a) Analytical model. The curve reaches an asymptotic value of D/d = 2. (b) Numerical model. The blue dots represent the slope at the origin of the mean square displacement ∆x 2 [see Fig. 2.13(b)]. The numerical results are in good agreement with the analytical model.

Figure 2 . 12 :

 212 Figure 2.12: Numerical model -(a) Sketch of the model. The grid is composed of a square lattice, and we show the example of a particle moving through the system (see text for the model parameters). (b) Results for 3 different paths (p z = 1/3). (c) Superposition of 4000 paths (p z = 1/3).A global morphology is obtained, where a 'diffusion' profile z ∼ x 2 is represented by the red line. Note that even if the representation of p z = 1/3 means that there is no hydrostatic pressure, the choice of having three directions only generates an asymmetric invasion profile and not a circular (isotropic) profile as expected.

2 >Figure 2 . 13 :

 2213 Figure 2.13: ∆x and ∆x 2 -(a)

Figure 2 . 14 :

 214 Figure 2.14: Cluster formation -The bubble position as function of time shows the stabilization around a constant value from time to time, revealing a cluster formation [h g = 14 cm, d =250-425 µm, Φ = 4.3 mL/s].

  . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Model description . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Morphology of air invasion: Bottom wall . . . . . . . . . .3.3.1 Geometrical characteristics of the invaded region . . . . . . . 3.3.2 Aspect ratio of the invasion zone . . . . . . . . . . . . . . . . 3.4 Air invasion in a free medium: Comparison . . . . . . . . 3.4.1 Extension and velocity of the invasion pattern . . . . . . . . 3.4.2 Curvature of the invasion zone . . . . . . . . . . . . . . . . . 3.5 Diffusion coefficient . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Mean square displacement ∆x 2 . . . . . . . . . . . . . . . . 3.5.2 Effective diffusion coefficient . . . . . . . . . . . . . . . . . . 3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ⊲ Morphology of air invasion in an immersed granular material G. Varas, V. Vidal & J.-C. Géminard, Phys. Rev. E 83, 061302 (2011).

Figure 3 . 1 :

 31 Figure 3.1: Images of the experimental pattern -For the same injected volume the path created by the injected gas has a drastically different geometry depending on the effective gravity g eff [V = 1.98 mL, (a) θ = 0 • , (b) θ = 80 • ].

Figure 3 . 2 :

 32 Figure 3.2: Model description and two examples of air paths from the numerical simulations -(a) Graphical sketch for the construction of the numerical model. We create a network which represents the capillary overpressure δP c that is distributed with random values taken from a normal distribution. We then add a vector that accounts for the hydrostatic δP g . Finally, the sum between these two terms gives the total pressure δP t at each vertex. (b,c) For the same injected volume (N = 4000), one observes that the geometry of the path depends drastically on the control parameter χ [(b) χ = 1000; (c) χ = 100]. As expected and also observed experimentally, for smaller χ (larger gravity effects) the air path exhibits less branches and reaches a larger height.

Figure 3 . 3 :

 33 Figure 3.3: Morphology of air invasion for the same injected volume V and different values of χ -The geometry of the invaded region goes from an almost circular pattern for a horizontal cell [(a) zero effective gravity, χ = ∞] to a vertically elongated shape when χ is decreased [(b)-(e)] (N = 8000, superposition of 2000 paths).

Figure 3 . 4 :

 34 Figure 3.4: Invasion density profile -(a) Surface plot of the density ρ p for the invasion at χ = 125. The results show an asymmetric density invasion in the vertical axis (z) with a maximum concentration at the injection point (z = 0 , x = 0) with ρ p = 1. (b) Alternative view for the density profile, the scale is limited to ρ p = 0.25 to visualize the outmost contour. Indeed, almost all the high densities are concentrated near the inlet [χ = 125, N = 8000].

Figure 3 . 5 :

 35 Figure 3.5: Analysis of the σ 2 -(a) Profiles of the probability density for three different altitudes z = 14, 65 and 240 [dashed line in Fig. 3.4 (b)]. The data can bed fitted by a normal distribution (black lines). (b) Square standard deviation σ 2 of the distribution as a function of the altitude z for different χ. The system exhibits a diffusive behavior for small χ. The curvature of the invaded region near the intlet is obtained from the slope at z = 0 (N = 8000). Inset: Zoom for small values of z and σ 2 and the initial slope for σ 2 vs z (red dashed lines).

Figure 3 . 6 :

 36 Figure 3.6: Shape of the invaded region vs. N and χ -Contours are reported in average over 2000 paths for f iso = 2 %. (a) Contours for increasing volumes (N = 2000, 4000, 6000 and 8000) for χ = 125. (b) Contours for decreasing effective gravity (χ = 1.6, 5.0, 12.5, 25, 50, 100, ∞) for the same injected volume (N = 2000).

Figure 3 . 7 :

 37 Figure 3.7: Contours, width ∆x and height ∆z of the invaded region and aspect ratio ∆z/∆x vs. χ -(a) We report the contours (lines) for f = 2, 4 and 6% (over 2000 paths, χ = 25 and N = 4000). (b) The width ∆x is defined as the width at half height, over the entire pattern of the intensity profile along the x-axis. (c) The height ∆z is that of the point at half the plateau value, over the entire pattern, of the intensity profile along the z-axis. (d) Aspect ratio ̺ ≡ ∆z/∆x vs. χ for three values of the injected volumes [ : N = 2000, : N = 4000 and • : N = 8000]. Inset: all the data collapse when reporting the aspect ratio as a function of χ/ √ N .

Figure 3 . 8 :

 38 Figure 3.8: Comparison of the simulations with and without bottom wall -(up) χ = ∞. (down) χ = 250. Even if the presence of a bottom wall changes one boundary condition, the global geometry is not altered by this restriction. (right) Representation for the ∆x/2 and ∆z/2 for both conditions. The red lines represent the results without bottom wall and the black lines with the bottom wall. For high χ the difference is about 10% while that for low χ the difference is of the order of 2%.

Figure 3 . 9 :

 39 Figure 3.9: Aspect ratio ∆z/∆x -Aspect ratio for the two types of bottom boundary conditions ̺ ≡ ∆z/∆x vs. χ/ √ N . The blue triangles represent the condition with bottom wall and the gray circle the condition without bottom wall. Note that all the data collapse on the same curve, following ρ ∼ χ -0.59 .

Figure 3 . 10 :Figure 3 . 11 :

 310311 Figure 3.10: Variation ∆z and ∆x as a function of the injected volume -(a) The typical vertical size of the structure ∆z increases linearly with N . (Inset) Normalized velocity as a function of χ. The velocity (the slope at the origin of the curves) ν z scales like d/ √ χ, the prefactor being of the order of the unity. (b) Typical horizontal size ∆x. Contrary to ∆z, in the horizontal direction the increase is slower and is an increasing function of χ. (Inset:) The normalized velocity ν x is almost constant in all the range of χ around a typical value of ν z /d ≃ 10 -2 .

Figure 3 . 12 :

 312 Figure 3.12: Curvature of the invasion zone -(a) Morphology of the invasion zone for χ = 250. (Zoom)The parameter l f is defined as the distance that separates the injection point (x i , z i ) with the bottom base of the resulting geometry. In red is represented the initial contour for f iso = 1% which is fitted by f (x) = x 2 /2a. (b) Ratio of l f /R c (0), with R c (0) = a (see text). The results show a transition at χ ∼ 10 3 . Over this value the invasion is close to a circular (non diffusive) invasion while under this value the parabolical shape accounts for the diffusive behavior.

  Figure 3.13: Normalized diffusion coefficient D/d vs χ -(a) Evolution of D/d obtained by considering the slope at the origin for the density profiles (see text). (b) Effective diffusion coefficient D eff /d. Both curves present a well defined slope until χ ∼ 10 3 .

Figure 4 . 1 :Figure 4 . 2 :

 4142 Figure 4.1: Final stationary state for different flow-rates Φ, and associated diffusion coefficient D -(a) Final invasion image produced by adding all the sequence of images [Φ = 0.91 mL/s]. The black lines represent the detected contour of the fluidized region. Note that the detection is limited until z ∼ 16 cm, thus avoiding possible deformation due to the free surface. (b) From the profile in Fig. 4.1 (a), we plot √ z vs x for both sides of the contour. The dashed line represents the slope of the curve, similar for both sides. (c) Results for all the flow-rates Φ (for one side of the invasion profile). The dashed line indicates the average slope. (d) Diffusion coefficient D vs Φ. The dashed line represents the average value of D = 2.12 ± 0.15 [h g = 22 cm, h w = 2 cm, d = 250 -425 µm, θ = 0 • ].

Figure 4 . 3 :

 43 Figure 4.3: Stable air-path structures -When tilting the cell for angles over 70 • the air path stabilizes and remains for long periods of time in the same configuration. Two examples for the same angle (θ = 80 • ) for a stable straight air-path (a) and a zig-zag air-path (b). At long times, the grains are pushed out of the granular bed, and it is impossible to define a fluidized zone.

Figure 4 . 4 :

 44 Figure 4.4: Evolution of the fluidized zone in time -(up) Generation of the fluidized zone. The deformation starts from the upper part of the granular bed (free surface) and in time, propagates toward the bottom. (down) Zoom over a region inside the fluidized zone. Bubbles are entrained down by the descending flow of grains. The speed of the flow decreases close to the border of the fluidized zone. The dashed arrows indicate the position of two bubbles along a period of time of one hour. Note that the bubble close to the center of the vertical of the injection point (the left bubble) descends faster [h g = 22 cm, h w = 24 cm, d = 250 -425 µm, φ = 0.91 mL/s, θ = 0 • ].

Figure 4 . 5 :

 45 Figure 4.5: Dynamics of the fluidized front -(a) Snapshot of the system. We analyze in a small region (zoom) the evolution of the fluidized front in the horizontal direction L f (x) and vertical direction L f (z) for different positions. (b) Evolution of L f (x) for different heights as a function of time. L f (x) is well defined after some minutes and then exhibits a logarithmic growth. Inset: Normalized L f (x)/ √ z. The data collapse when normalized by the square root of the height z. (c) Evolution of L f (z) for different locations x. The data share the same asymptotic slope.

Figure 4 . 6 :

 46 Figure 4.6: Evolution of the stick-slip zone -(a) At the beginning of the injection we observe three different zones. The first is the region where the grains flow (fluidized zone).Next to this, the stick-slip zone (SSZ), is characterized by spontaneous events. Finally, the grains in the outer region are compact and motionless. (b) In time, the fluidized zone invades the stick-slip zone which finally disappears. Note that the compact zone is invariant along all the experiment (in some cases the only alteration of this zone is observed during the initial injection, especially for the percolation regime) [h g = 22 cm, h w = 2 cm, d =250-425 µm, Φ = 2.17 mL/s].

  4.7). (a) Space-time diagram for the surface height evolution (zero is the initial flat surface). (b) Space-time diagram of the free surface amplitude variation obtained by subtracting two consecutive images. (c) Binarization of (b). In each figure, the color scale indicates the variation in height expressed in centimeters.

Figure 4 . 7 :

 47 Figure 4.7: Space-time diagram for the evolution of the free surface of grains -

Figure 5 . 1 :

 51 Figure 5.1: Crater formation -Snapshot of the system (the image color is inverted). A series of ascending rolls pushes the grains away from the vertical central axis. Subsequently, the grains gently deposit back onto the free surface of the bed. Along the piles flanks, provided that the local angle exceeds the angle of avalanche, we observe continuous avalanches. Inside the crater, the flowing granular material partly replaces the grains advected upwards at the center. Along the outer flanks, the deposited granular material either flows or sits at the free surface. Inset: Zoom for the outer side of the left pile. The dashed line is a guide for the eye that indicates the separation between the initial free surface of grains and the new material deposited [d = 400 µm, h g = 8 cm, h w = 16 cm, Φ = 3.5 mL/s].

Figure 5 . 2 :

 52 Figure 5.2: Overview of the gas flow process -(a) At intermediate flow-rate Φ, the pressure signal exhibits a spontaneous alternation between activity and rest periods which correspond respectively to the bubbling and open-channel regimes. The emission of the bubbles at the free surface accordingly exhibits two qualitatively different behaviors : (b)In the bubbling regime, a large bubble forms underneath the free surface of the granular bed and then bursts. The explosion pushes out of the bed a large number of grains, which are subsequently advected upwards in the bubble wake. (c) In the open-channel regime, the walls of the channel are stable and a continuous gas flow escapes the granular bed. In this regime the grains are only torn out from the granular bed by the water flow and subsequently advected in the wake of the small bubbles that form in water [d = 400 µm, h g = 8 cm, h w = 16 cm, and Φ = 2.5 mL/s].

Figure 5 . 3 :Figure 5 . 4 :

 5354 Figure 5.3: Temporal evolution of the crater -We define L, the distance between the two peaks, and x, the distance from the center. The angle α and β denote, respectively, the maximum angle that the inner and outer flanks make with the horizontal [d = 400 µm, h g = 8 cm, h w = 16 cm, and Φ = 3.5 mL/s].

Figure 5 . 5 :

 55 Figure 5.5: Inner and outer angles evolution -(a) Definition of α and β. The red circles denote the vertical to measure the angle. Note that for large x, the slope of the outer flank changes. (b)The grey zone marks the early stages of the crater growth during which the crater size compares with the size of the bubbles. The angle α, which is then not measurable, remains constant during the whole experimental time. In contrast, the angle β increases continuously to reach a constant value after a finite time. In order to show that the slope is not constant along the outer flanks, we report the temporal evolution of the slope at x = 1 2 L(t) + 1.5 cm: the local slope remains smaller than β until reached by the avalanche flow. Lines are only guides for the eye [Φ = 1.67 mL/s, d = 400 µm, h g = 8 cm, h w = 16 cm].

Figure 5 . 6 :

 56 Figure 5.6: Effect of the water height h w on the crater growth -(a) Distance L as a function of time for different water heights h w . The dashed lines shows the distance L for t = 30 min. (b) Distance L(30 min) vs water height h w . For h w < L/4, the growth of the crater is limited by the water height and L(30 min) increases linearly with h w . In contrast, for h w > L/4 (blue region) the growth of the crater is not altered by the water height and L(30 min) does not depend on h w [d = 100 µm, h g = 8 cm, Φ = 3.06 mL/s].

Figure 5 . 7 :

 57 Figure 5.7: Rolls size -(a) Snapshot of the system (the image color is inverted) where the zone delimited by the rolls is represented. The distance L r is measured from the center to the border, hence, the total distance reached by the rolls will be 2 × L r . (b) Comparison of L r with the crater size L. The free surface, and hence the crater, does not affect the size of the rolls L r . The results shows that there is no correlation and L r is constant in time [d = 400 µm, h g = 8 cm, Φ = 3.06 mL/s].

5. 4 Figure 5 . 8 :

 458 Figure 5.8: Distance L vs. time t -The distance L between the summits increases

Figure 5 . 9 :

 59 Figure 5.9: Grain size dependence -Snapshot of the system (the image color is inverted) for two different grain size [(up) d = 100 µm, (down) d = 400 µm].The volume of grains lifted by the bubbles is drastically different from one grain size to the other. For the small grain size (up) there is almost a suspension that contain larger quantities of grains, while for the large grain size (down) there are two welldefined convective rolls with less grains in suspension. As expected the growth of the crater for the same flow-rate Φ is faster for the small grain size [Φ = mL/s].

Figure 5 . 10 :

 510 Figure 5.10: Height evolution in time -(a) Local height h(x, t) vs time t. Far away from the center, the height h of the granular layer increases due to the deposition of the grains. We observe a linear increase of h(x, t) with time until the region is reached by the grains avalanching along the outer flank. (b) Local growth velocity ∂h/∂t(x, t) vs distance x. The data are successfully interpolated by a x Lc exp (-x Lc ) for large x, and we get at a rough estimate L c = (1.13 ± 0.04) cm [Φ = 2.93 mL/s, h g = 16 cm, h w = 8 cm].

Figure 5 . 12 :Figure 5 . 13 :

 512513 Figure 5.12: Rayleigh-Taylor instability -After tilting the cell, the granular bed is above the water layer. Then, the instability is triggered and evolves until reaching the top surface. Granular fingering, mixing and sedimentation are observed.

Figure A. 1 :

 1 Figure A.1: Experimental setup -A film (dashed line) is drawn at the open end of a cylindrical tube (black lines).Air is injected in the system with a syringe and the bubble (blue line) grows. The arrows indicate the external forces applied to the system (pressure and surface tension). A zoom on the film shows the bubble interface. We consider a constant thickness e, R is the radius of curvature, r the cylinder radius, V the volume of the bubble, P 0 the atmospheric pressure and ∆P the overpressure inside the tube.

Figure A. 2 :

 2 Figure A.2: Interference fringes and film thickness evolution -(a) Interference fringes showing five peaks. The number of peaks and the difference between them are used to calculate the thickness (Eq. A.1). (b) Film thickness evolution until its rupture. Note that it is difficult to quantitatively reproduce these results (see text).

Figure A. 3 :

 3 Figure A.3: Film bulging when injecting a constant air flow-rate -The film geometry is later fitted by an circle arc (see text) [Φ = 25 mL/s].

Figure

  Figure A.4: Comparison between the theoretical and experimental overpressure

Figure A. 5 :

 5 Figure A.5: Pressure evolution and radius of curvature of a bubble -(a) Pressure evolution for the different flow-rates as a function of the normalized injected volume Φ t/V 0 [Φ = 0.5 -26 mL/s]. The results show a similar behavior for each curve but a different maximum pressure. Note that the initial slope is the same for all the data, and the separation occurs near and beyond the maximum. Inset: Zoom on the region near the maximum. (b) Maximum pressure as a function of the flow-rate. The data reveal that there is a linear correlation between the maximal pressure and the flow-rate. (c) Radius of curvature as a function of time for different flow-rates [ = 6 mL/min, • =10 mL/min , △ = 20 mL/min and ⋄ = 26 mL/min]. All the curves display the same minimum R = 8.25 mm which corresponds to the external radius of the tube. (d) Overpressure ∆P as a function of 1/R. The slope of the curve is the surface tension γ.

Figure

  Figure A.5 (b) shows the maximum pressure value ∆P max for each flow-rate. The results were calculated by taking the average of three series, with the error bar representing the standard deviation. The data reveal that there is a linear correlation between the maximal pressure and the flow-rate which could be the consequence of a dynamic surface tension. Note that the first fourth data points were taken with the same sample but another day. This could be the at the origin of the discrepancy with the global linear trend.

FigFigure A. 6 :

 6 Fig. A.5 (c) displays the radius of curvature as a function of time for different flow-rate Φ. As expected all the curves have the same minimum value corresponding to the radius of the tube and the maximal pressure. The increase after this point depends on the air flow. In Fig. A.5 (d) we report the overpressure as a function of the inverse of the radius of curvature 1/R. The slope of the curve gives the surface tension. The

  1 mL/ s͒ by means of a home-made flow meter. A sensor ͑MKS Instruments, 223 BD-00010 AB͒, connected to an acquisition board ͑National PHYSICAL REVIEW E 79, 021301 ͑2009͒ 1539-3755/2009/79͑2͒/021301͑7͒ ©2009 The American Physical Society 021301-1

  FIG. 4. Temporal evolution of the crater geometry: ͑a͒ Inner and outer angles ␣ and ␤ vs time t ͑⌽ = 1.67 mL/ s͒. The grey zone marks the early stages of the crater growth during which the crater size compares with the size of the bubbles. The angle ␣, which is then not measurable, remains constant during the whole experimental time. In contrast, the angle ␤ increases continuously to reach a constant value after a finite time. In order to show that the slope is not constant along the outer flanks, we report the temporal evolution of the slope at x =

1 2

 1 FIG.5. The granular flows: The second configuration is used: We light up the sample from top with the linear light source and take pictures with the high-resolution digital camera. Initially, a thin layer of colored grains is deposited at the free surface of the initially flat and horizontal bed. Then, a series of ascending rolls pushes the grains away from the vertical central axis. Subsequently, the grains gently deposit back onto the free surface of the bed. Along the piles flanks, provided that the local angle exceeds the angle of avalanche, we observe continuous avalanches. Inside the crater, the flowing granular-material partly replaces the grains advected upwards at the center. Along the outer flanks, the deposited granular material either flows or sits at the free surface. One can clearly observe, in the inset, that new material, deposited far away from the center forms a thin layer of material which remains at rest. The deposition flow rate, Q͑x͒, is estimated from the temporal evolution of the growth velocity ץh / ץt of the deposited layer in that region ͓Fig. 4͑c͔͒ ͑d = 400 m, h g = 8 cm, h w = 16 cm, and ⌽ = 3.5 mL/ s͒.

B. 1 .

 1 Influence of the gas flow ⌽ and grain size d Dependence on the air flux ⌽

FIG. 6 .

 6 FIG.6. Distance L͑30 min͒ vs water height h w : For h w Ͻ L / 4, the growth of the crater is limited by the water height and L ͑30 min͒ increases linearly with h w . In contrast, for h w Ͼ L / 4, the growth of the crater is not altered by the water height and L ͑30 min͒ does not depend on h w ͑d = 100 m, h g = 8 cm, and ⌽ = 3.06 mL/ s͒.

FIG. 7 .

 7 FIG. 7. Distance: L vs time t:The distance L between the summits increases logarithmically with time, independently from the gas emission regimes. Indeed, at small ⌽, the gas flow regime consists mainly in the emission of independent bubbles, whereas it consists mainly in a continuous emission of gas at large ⌽. From the interpolation of the experimental data with L / L 0 =ln͑⌽t / V 0 ͒, we estimate L 0 = ͑3.5Ϯ 0.7͒ cm and V 0 = ͑0.4Ϯ 0.2͒ mL for d = 400 m ͑black curves͒ whereas we estimate L 0 = ͑3.5Ϯ 1͒ cm and V 0 = ͑0.02Ϯ 0.01͒ mL for d = 100 m ͑gray curve͒ ͑h g = 8 cm, h w =16 cm͒.

FIG. 1 .

 1 FIG. 1. Experimental setup:Air is injected at a constant flow rate into a chamber of volume V connected to the bottom of an immersed granular column (h g is the height of the granular column and h w is the height of the water above). The loci of the gas emission are determined by imaging the free surface by means of a webcam.

FIG. 2 .

 2 FIG. 2. Loci of the gas emissions: Each of the white circles indicates the locus of a bubble emission. Half of the bubbles were emitted in the black circle of radius ζ (Fig. 3). Inset: Example of a bubble emission (d =250-425 μm, = 3 mL/s, h g = 20 cm).

r

  0 p(u)du[Fig. 3(c)]. We observe that F is rather well described by a Gaussian function [Fig. 3(c), inset]. We define the typical width ζ of the central region by F (ζ ) ≡ 1/2, which means that half of the bubbles escape the free surface in the circle of radius ζ at the vertical of the inlet.

  FIG. 3. (a) Distance r vs time t (data from Fig. 2). A total acquisition lasts for 10 min (300 s are displayed here). (b) Autocorrelation R of the signal as a function of the neighboring bubble number n. (c) Complementary cumulative distribution function 1 -F (r). The distribution is described by a Gaussian function (full line). The width at half height, 2ζ , corresponds to the diameter of the circle in which half of the bubbles are emitted. Inset: Log-lin plot of the data (circles) and Gaussian fit (full line) (d = 250-425 μm, = 3 mL/s, h g = 20 cm).

FIG. 5 .

 5 FIG. 5. Air paths observed in the 2D experiment. The air creates several paths within the granular bed and one observes several gasemission loci at the free surface (arrows) (h g = 22 cm, h w = 2 cm, d = 250-425 μm, = 0.43 mL/s).

FIG. 6 .

 6 FIG.6. Image of the fluidized region. The contrast results from the lower grain compaction in the fluidized region. Its shape is determined by image analysis, and subsequently interpolated by a parabolic function (white dashed line). We denote by 2ζ * (z) the characteristic width, in the bulk, of the region explored by the air in the stationary regime at height z (h g = 22 cm, h w = 2 cm, d = 250-425 μm, = 2.17 mL/s).

  FIG. 7. (Color online) (a) Definition of the capillary pressures δP associated with the possible directions of air path after reaching a node (black dot): toward x < 0 (δP -), x > 0 (δP + ), or z > 0 (δP z ). The model discards the possibility for the air to propagate downward. (b) Sketch of the considered grain arrangement. The air, injected at the bottom, propagates through the easiest path between the grains (dashed line), biased by gravity (see text). (c) Normalized coefficient D/d vs parameter χ . The coefficient D scales like the grain size and is a function of the unique parameter χ [Eqs. (1) and (4)].

  )Thus, a diffusionlike process accounts for the dependence of ζ * on z and that of ζ on h g . The associated coefficient D scales like the grain size d. However, to account for the dependence 011302-4 of D on d, one must evaluate the dependence of the probability p z on the problem parameters.

FIG. 8 .

 8 FIG.8. Radius r vs bubble number n. For large depth, a strong correlation between successive bubbles might be observed: several successive bubbles are emitted in a narrow region (cluster) before the system turns back to the regular regime (h g = 14 cm, d = 250-425 μm, = 4.3 mL/s.)

  FIG. 1. (Color online) Sketch of the experimental setup. A total volume V of air is injected at constant flow rate into a thin vertical cell (Hele-Shaw cell) containing a granular layer immersed in water.The effective gravity is changed by tilting the cell by an angle α (from 0 • to 90 • ). The formation of the pattern in the granular bed is recorded with a camera positioned in front of the experimental cell.

FIG. 2 .

 2 FIG.2. Images of the experimental pattern. The 2D experiment reveals that, for the same volume V = 1.98 mL, the path created by the injected gas has a drastically different geometry depending on the effective gravity g eff . (a) For an almost horizontal cell (α = 80 • ), one observes a highly branched pattern, which remains located around the outlet (note that its typical width compares with its typical height). (b) By contrast, for a vertical cell (α = 0 • ), even if the image reveals a significant number of branches, the maximum height reached by the gaseous structure is much larger than its width (d = 250-425 μm).

FIG. 3 .

 3 FIG. 3. Air paths from the numerical simulations. Note first that the numerical procedure allows the creation of side branches. For the same injected volume, thus the same number N = 4000 of occupied vertices, one observes that the geometry of the path depends drastically on the control parameter χ [(a) χ = 1000; (b) χ = 100].As expected and also observed experimentally, for smaller χ (larger gravity effects) the finger exhibits less branches and reaches a larger height.

FIG. 4 .

 4 FIG. 4. Superposition of 2000 paths for the same injected volume V and different values of χ . The morphology of the invasion region goes from an almost circular pattern for a horizontal cell [(a) zero effective gravity, χ = ∞] to a vertically elongated shape when χ is decreased [N = 8000, (a) χ = ∞; (b) χ = 1250; (c) χ = 250; (d) χ = 166; (e) χ = 100].

2 FIG. 5 .

 25 FIG. 5. (Color online) Contours, width x, and height z of the invaded region and aspect ratio z/ x vs. χ . (a) We report the contours (lines) for f = 2,4, and 6 ‰ (over 2000 paths, χ = 25 and N = 4000). (b) The width x is defined to be the width at half height, over the entire pattern, of the intensity profile along the x axis. (c) The height z is that of the point at half the plateau value, over the entire pattern, of the intensity profile along the z axis. (d) Aspect ratio ≡ z/ x vs. χ for three values of the injected volumes ( : N = 2000, : N = 4000 and • : N = 8000). Inset: A collapse of the measurements is observed when reporting the aspect ratio as a function of χ/ √ N .

FIG. 6 .

 6 FIG. 6. Shape of the invaded region vs. N and χ . Contours are reported on average over 2000 paths for f = 2 ‰. (a) Contours for increasing volumes (N = 2000,4000,6000, and 8000) for χ = 125. (b) Contours for decreasing effective gravity (χ = 1.6,5.0,12.5,25,50,100,∞) for the same injected volume (N = 2000).

FIG. 7 .

 7 FIG. 7. Analysis of the probability profiles. (a) Superposition of 2000 paths for χ = 166. The gray horizontal lines are associated to the profiles reported in the right panel. The dashed gray lines represent three different heights. (b) Profiles of the probability density at three different altitudes z = 14,65, and 240. (c) Square σ 2 of the distribution width as a function of the altitude z for χ = 5,50,166,500,1666, and ∞. The system exhibits a diffusive behavior for small χ . The curvature of the invaded region near the outlet is obtained from the slope at z = 0 (N = 8000).

FIG. 8 .

 8 FIG.8. Coefficient D eff and velocity ν z vs. parameter χ . The radius of curvature of the invaded region, near the outlet, is independent of the injected volume ( : N = 2000, : N = 4000, and • : N = 8000) and is given by D eff d χ provided that the gravity effects are large enough (χ 30 √ N ). In the same conditions, the typical vertical size of the structure z increases linearly with N . The velocity (the slope) ν z scales like d/ √ χ , the prefactor being of the order of unity.

Table 1 . 1 :

 11 Grain size distribution analysis for the different batchs -The left column indicates the batch (from USF Matrasur).

	grain size d (µm) mean diameter µ (µm) standard deviation σ (µm)
	200 -250 250 -425 400 -500	203.3 ± 0.8 293 ± 3	19.3 ± 0.4 50 ± 6

Table 2 .

 2

	• C)	ρ a	1.204	kg/m 3
	Permeability	k	100	mD

1: Parameters in the experiments and their associated dimension

Table A .

 A 1: Typical values for a soap film at room temperature.

	A.3 Dynamic pressure evolution	123
	Parameter Value Units
	γ	10 -2
		8)

  Germán Varas, Valérie Vidal, and Jean-Christophe Géminard Laboratoire de Physique, Université de Lyon, Ecole Normale Supérieure, CNRS, 46 Allée d'Italie, 69364 Lyon Cedex, France ͑Received 6 November 2008; published 2 February 2009͒

  In 3D (triangles), the typical size ζ of the region in which bubbles escape the free surface is compatible with the scaling ζ ∝ h g , independent of the flow rate . In 2D (circles), accordingly, ζ * scales like

	√	z, almost
	independent of [gray scale, (mL/s): white, 1.0; gray, 2.0; black,
	3.0 (3D experiment); white, 0.4; gray, 0.9; black, 1.3 (2D experiment);
	d = 250-425 μm.]

1/2 ) FIG. 4. Distances ζ vs h g and ζ * vs √ z.

  Germán Varas, Valérie Vidal, and Jean-Christophe Géminard Laboratoire de Physique, Université de Lyon, Ecole Normale Supérieure -CNRS, 46 Allée d'Italie, F-69364 Lyon Cedex, France (Received 19 April 2011; published 23 June 2011)

  ).

	DOI: 10.1103/PhysRevE.83.061302	PACS number(s): 83.80.Fg, 47.57.Gc, 47.85.Dh
	I. INTRODUCTION

The Toarcian is a period or stage in the Early or Lower Jurassic. It spans between 183.0 Ma (million years ago) and 175.6 Ma.

Although a more specific classification can be made for imbibition, for example: spontaneous imbibition, constant influx imbibition, quasistatic imbibition and dynamic invasion with constant flow-rate of the displacement fluid[START_REF] Payatakes | Immiscible microdisplacement and ganglion dynamics in porous media[END_REF][START_REF] Sahimi | Flow phenomena in rocks -from continuum models to fractals, percolation, cellular-automata, and simulated annealing[END_REF].

The hydraulic gradient, called also Darcy slope, is defined by the difference of height of water ∆h over a length l.

A fourth regime was reported by the authors, but it corresponds to a transition regime and is not described here.

A Hele-Shaw cell is a quasi-two-dimensional cell where the width is large compared to the cell thickness.

One can distinguish between colloidal suspensions where thermal agitation of the grains is important and granular suspensions where thermal agitation is negligible.

The phenomenon was observed by[START_REF] Zoueshtiagh | Effect of a vertically flowing water jet underneath a granular bed[END_REF] but not studied in their work.

We only present here three batchs characterization. Different grains size were used in this thesis and all the batchs are characterized by a Gaussian distribution.

All the parameters and the range used are summarized in table

2.1.1 of this chapter.

As a consequence of the grain volume conservation, small displacements of the free surface are observed.

We only interpolate the bulk region of the fluidized zone, excluding the upper part that, as a consequence of the grain-water interface, exhibits a deformation close to the free surface.

polydispersity and packing fraction.

In order to change any of these parameters we would have to change the medium completely which is long and tediously. For more information see Sec.

2.1.2.

In the experiment, even if we control carefully the different parameters (d, Φ, ...), the local packing can vary when we prepare the initial condition (see Sec. 2.1.2) and the air path will thus differ.

For some simulations, especially for lower χ the vertical size will be higher in order not to let the air path reach the surface.

The profiles are well fitted by a normal distribution until z <

Above this value (specially for greater χ), the lack of sufficient data limits the analysis

We will discuss this point in the next section.

This is a perfect circle, considering that the statistics is restricted by the limited number of repetitions.

Morphology of air invasion: Numerical analysis

As the compacted zone remains fixed in time, the superposition of all the images generates a large contrast difference in the final image.

For sake of clarity, we only displayed the points associated with one side for the different profiles.

The reason of this dark color is unknown to us. Maybe it is due to partial segregation and/or some dust which is trapped inside the cell.

This may not be confused with the cluster formation shown in section 2.6. In a quasi-2D experiment it is more easy for a path to stabilize, resulting in a fixed bubble emission.

 5 Note that the results could be improved with a better temporal and spatial resolution.

It is for this reason that we added all the images to find the final state of the fluidized zone (Sec. 4.1).

The correlation between successive bubble emission has been analyzed in Sec.

2.2.2. Similar results are found here.

The existence of the bubbling and open-channel regimes and the transition between them are tightly related to the non-Newtonian rheology of the granular material as proven by an experimental study conducted in another non-Newtonian material, exhibiting the same rheological characteristics (existing of a yield stress and shear thinning)[START_REF] Divoux | Acoustic emission associated with the bursting of a gas bubble at the free surface of a non-Newtonian fluid[END_REF].

The outer flank exhibits a bump shape.

The precise dependence of L 0 and V 0 on d is difficult to access experimentally as such measurements require a precise knowledge of the origin of time (a delay alters significantly the slope and the offset in the semilog plot) which we are missing.

As the model is in two dimensions, v g is in fact an area. We still maintain the notation of a volume to be consistent with the experiment.

Conclusion & Perspectives

The spectrometer measures the light reflected by the film, which varies according to its thickness. Note that the film thickness estimation is limited to the region analyzed by the spectrometer (top of the bubble).

This measurement requires further improvement. For example the vertical position of the bubble changes in time while the spectrometer remains fixed at the same position. This effect restrains the analysis for large gas volume injected.

Light-Emitting Diode.

As already presented above this hypothesis is valid as the bubble (at least at the top) maintains an almost constant thickness for several seconds before it ruptures.

It is difficult experimentally to vary Φ over more than a decade.

Remerciements

ACKNOWLEDGMENTS G.V. acknowledges a grant by CONICYT (Comisión Nacional de Investigación Científica y Tecnológica, Gobierno de Chile). The authors are grateful to F. Vittoz, from ENS Mechanical Workshop, for the design and setup of the 2D experiment. ACKNOWLEDGMENTS G.V. acknowledges a grant by CONICYT (Comisión Nacional de Investigación Científica y Tecnológica, Gobierno de Chile). The authors thank F. Vittoz for the technical support.

Appendix A Dynamic surface tension of bursting bubbles A.1 Introduction

Understanding the physics of bubble bursting is of importance in fundamental physics as well as in natural phenomena. The bursting of bubbles has been studied for different types of liquid, for example soap films, viscous liquids, crystal liquids, etc. [START_REF] Bergeron | Forces and structure in thin liquid soap films[END_REF][START_REF] Debregeas | The life and death of "bare" viscous bubbles[END_REF][START_REF] Mueller | Experimental study of the bursting of inviscid bubbles[END_REF]. Other interesting phenomena have been observed, for example, the vortices produced by thermal convection at the surface of a bubble that resemble hurricanes [START_REF] Seychelles | Thermal convection and emergence of isolated vortices in soap bubbles[END_REF]. In the geophysical context they play an important role in volcanic eruption, one of the main objectives in this area being the understanding of the explosion of large gas bubbles at the top of the volcanic conduit. During the last decades monitoring has been performed on several volcanoes where large gas bubbles rise and explode at the top of the magma conduit [e.g. [START_REF] Jaupart | The generation and collapse of a foam layer at the roof of a basaltic magma chamber[END_REF]; [START_REF] Vergniolle | Dynamics of degassing at Kilauea Volcano, Hawaii[END_REF]; [START_REF] Vergniolle | Origin of the sound generated by strombolian explosions[END_REF]]. These explosions generate both seismic and acoustic waves propagating in the ground and in the atmosphere, respectively. The huge overpressure inside the bubble before bursting (∼ atm), inferred from numerical models [Vergniolle and Brandeis (1996); Vergniolle et al. (1996)], is responsible for the violence of the explosion. To comprehend this phenomenon, experiments aiming at characterizing the acoustic emission associated with bursting bubbles have been done in viscous liquids [START_REF] Jaupart | Laboratory models of hawaiian and strombolian eruptions[END_REF] as well as in non-Newtonian fluid [START_REF] Divoux | Acoustic emission associated with the bursting of a gas bubble at the free surface of a non-Newtonian fluid[END_REF]. However, no model yet is able to explain the origin of such overpressure (usually MPa).

In this appendix, motivated by the will to understand of these natural phenomena we report preliminary work on the study of bubble dynamics. The idea behind this work is to find if any correlation exists between the dynamics of the bursting bubble and the pressure at bursting. In order to understand what determines the overpressure inside a bubble at bursting, and the role of the fluid rheology, we performed laboratory experiments. A bubble of well-controlled geometry is mimicked by a thin liquid film [2]. Nous présenterons quelques résultats préliminaires montrant la dynamique d'évolution de la zone fluidifiée en fonction des différents paramètres. Enfin, lorsqu'on augmente la hauteur d'eau au-dessus de la couche granulaire, l'advection et le dépôt successif des grains forment un cratère, composé de deux dunes qui croissent et s'éloignent du centre. La taille typique du cratère augmente de façon logarithmique dans le temps, indépendamment du processus d'émission du gaz [3].

Air rise through an immersed granular bed -bulk and surface dynamics

When air is injected at the bottom of an immersed granular layer, it crosses the system by percolating or fracturing. It thus forms several paths that reach the free surface of the layer at different locations. In this thesis, we study this process experimentally (for a three and two dimensional setup), numerically and theoretically. First, we focus on the dynamics of the air invading the medium at short and long time scale, when injecting a continuous air flow. At long time, the typical size of the region explored by the air can be accounted for by a diffusion-like process [1]. We also investigate the effect of gravity by tilting the experimental cell. We contrast the results with numerical simulations for the injection of a fixed volume of air, and characterize the morphology of the invasion zone. We show that the typical height and width of the region explored by the air does not depend on the injected volume only, but also on a dimensionless parameter χ which accounts for the relative effects of the gravity and capillarity [2]. Finally, when increasing the water height above the granular layer, successive grain advection and deposition form a crater consisting of two dunes growing and moving apart one from the other. We observe that the typical size of the crater increases logarithmically with time, independently of the gas emission process [3].