# Unifying fluidization and defluidization of granular columns

L Chupin, Olivier Roche, Siet van den Wildenberg

► To cite this version:

L Chupin, Olivier Roche, Siet van den Wildenberg. Unifying fluidization and defluidization of granular columns. Biophysical Reviews, 2025, 111 (4), pp.045421. 10.1103/PhysRevE.111.045421. hal-05043503

HAL Id: hal-05043503

https://uca.hal.science/hal-05043503v1

Submitted on 23 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Unifying fluidization and defluidization of granular columns

L. Chupin,¹ O.Roche,² and S. van den Wildenberg²,³,*

¹Université Clermont Auvergne, CNRS UMR 6620,
Laboratoire Mathématiques Blaise Pascal, F-63000 Clermont-Ferrand, France

²Université Clermont Auvergne, CNRS, IRD, OPGC,
Laboratoire Magmas et Volcans, F-63000 Clermont-Ferrand, France

³Université Clermont Auvergne, CNRS, Laboratoire de Physique de Clermont Auvergne, F-63000 Clermont-Ferrand, France

(Dated: April 22, 2025)

We study air pressure during the fluidization and defluidization of granular columns with the aim to unify the description of these processes. To this end, we perform experiments with different air flow velocities and columns of different heights. The experimental measurements evidence a non-linear relationship between pressure and the depth in the column. This complex pressure profile challenges the conventional framework of Darcy's law, typically used to analyze airflow in granular columns. We introduce a novel law that effectively captures these non-linear effects, occurring even at low particle Reynolds numbers. An advantage of our approach is its applicability to defluidization—a key phenomenon in gravity-driven gas-particle flows. Using this new law, we present a model to describe the decrease of air pressure in a defluidizing column, with numerical simulations aligning well with experimental results.

## INTRODUCTION

The interstitial gas in granular mixtures plays an important role in many industrial and geophysical processes. If the pore gas pressure is high, the effective friction is reduced and the flow of such mixtures is fluid-like. Instead, if the pore pressure is negligible the mixture will flow like a dry granular material [1]. To accurately predict the behavior of gas-particle mixtures, it is thus important to understand the gas solid interactions. The basic law governing the laminar flow of a fluid through a porous media, such as an ensemble of grains, is Darcy's law. If the flow is directed in the direction denoted z then this law predicts a proportionality between the pressure gradient $\partial_z P$ and the velocity U of the flow in this direction:

$$ \frac{\partial P}{\partial z} = -\frac{\mu}{k} U. \qquad (1) $$

The coefficient depends on the dynamic viscosity of the fluid $\mu$ and the permeability of the granular medium $k$. Note that U is the superficial fluid velocity, or, in other words, a macroscopic averaged velocity. For flows through porous media with higher particle Reynolds numbers, inertial effects may become significant. In this flow regime, the Darcy-Förchheimer equation captures the nonlinear behavior of the pressure gradient as a function of flow velocity [2, 3]:

$$ \frac{\partial P}{\partial z} = -\frac{\mu}{k} U + \frac{\rho}{k_1} U^2 \qquad (2) $$

with the fluid density $\rho$ and the inertial permeability $k_1$. For monosize beads, the permeabilities in the viscous and inertial regime can be obtained from the solid volume fraction $\phi$, and the diameter $d$ and sphericity $\theta$ of the particles, see [4]: $k = \frac{\theta^2 d^2 (1 - \phi)^3}{150 \phi^2}$ and $k_1 = \frac{\theta d (1 - \phi)^3}{1.75 \rho \phi}$.

The fluidization of granular media is principally characterized via the minimum fluidization velocity $U_{mf}$, which is defined as the superficial fluid velocity at which there is a transition from a fixed bed to a fully fluidized state. Multiple studies have developed empirical and semi-empirical correlations of $U_{mf}$ [5–9].

In contrast to fluidization, the defluidization process, in which the fluid escapes from the granular bed, is much less studied. Nevertheless, defluidization is important in many natural systems since the fluid is not continuously supplied. An important question is thus how the fluid pressure in such systems decreases over time. Indeed, it was shown that diffusion laws can quantitatively describe the time evolution of pore fluid pressure, but the inferred diffusion coefficient increases with the height of the granular column [1, 10]. It was suggested that this effect is simply due to the extra volume of air below the bed in the so-called wind-box [10]. However, as we will show, it occurs even when the extra volume is negligible compared to the volume of interstitial gas in the bed.

We performed careful experiments on the fluidization and defluidization of granular columns of different heights. We show that the fluid pressure profile increases nonlinearly inside the granular bed contrary to predictions by Darcy's law. We capture this effect by extending Darcy's law with a linear pressure term. This simple extension permits us to easily rewrite the law in terms of the fluid flow and test our model in defluidization experiments. We find good agreement between experimentally measured pressure decays and the extended Darcy's law.

## EXPERIMENTAL METHODS

Experiments were performed in a plexiglass cylindrical pipe with inner diameter $D = 0.14$ m and a height of 2 m
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Schematic diagram of the experimental set-up, a cylinder placed on a the porous plate W and a windbox A. Pressure sensors $S_i$ are located at positions $z_1 = 0.04$ m, $z_2 = 0.13$ m and $z_3 = 0.3$ m. Air flow is controlled by a solenoid valve V.(b) Particle size distribution with $d$ the equivalent circle diameter (diameter of a circle with the same area as the projected area of the particle image) and $P_v$ the volume fraction in each bin. (c) Solid volume fraction $\phi$ versus the height of the column $H$.</figcaption></figure>

(Fig. 1a). At the bottom of the pipe was a porous plate (W) and the ensemble was placed on top of a chamber called the wind-box (A) with the same inner diameter as the pipe. The porous plate had a thickness of 5 mm and a pore size of 20 µm. The wind-box and the porous plate are essential to ensure uniform airflow within the granular column. The "extra volume" of air refers to the air trapped in the wind-box (between the valve and the granular column). When the valve is closed, this trapped air must escape through the column, which can delay and alter the pressure evolution if its volume is significant compared to the air volume in the column. Therefore, it is important to minimize the wind-box volume, which we achieved by filling it with glass beads ($d = 1.5$ mm). The volume of air in the wind-box and the porous plate was estimated to be approximately $4.8 \times 10^{-5}$ m³, which corresponded to about 0.1% of the volume of air in the smallest granular column tested here. Taking a permeability of the porous plate of $k_{plate} \approx 10^{-11}$ m² and $\phi_{plate} \approx 0.5$, the pressure drop over the plate was estimated to be of order of 52 Pa for the highest air flow tested.

A laboratory air supply system delivered mean air flow velocities (U) through the wind-box and the porous plate into the pipe. The pressure of the interstitial air in the cylinder was measured by three piezo resistive sensors $S_1, S_2$ and $S_3$ installed at a distance from the porous plate of 0.04 m, 0.13 m and 0.3 m, respectively. The sensors were placed in adapters so that the front of the sensors was in direct contact with interstitial air via a small

hole covered by a grid (36 µm mesh size) to prevent particles from entering the adapter and colliding with the sensors. The back of the sensors were in contact with the atmospheric pressure only. Therefore, the sensors measured the differential pressure between their front and back sides, which corresponded to the pore fluid pressure in the pipe.

We used nearly spherical glass beads of density $\rho = 2500 \text{ kg.m}^{-3}$, which were characterized using MORPHOLOGY G3SE$^{TM}$ equipment from Malvern Instrument$^{TM}$ following the protocol of Leibrandt and Le Pennec [11]. The particle size distribution has a tail at smaller particle sizes and, therefore, we consider the Sauter diameter. The Sauter diameter accounts for the importance of drag forces on the particle [12]. For the glass particles, we find a Sauter diameter $d = 62$ µm and a sphericity $\theta = 0.97$.

To generate reproducible granular packings we developed a sample preparation protocol as follows. First, a known mass of glass beads was poured into the pipe. Second, to avoid excessive cohesion due to ambient moisture, the glass beads were aerated for about 20 min using a dry air flux at a velocity greater than $U_{mf}$. Third, the granular bed was expanded using the air flux, after which the air flux was stopped and the bed was compacted by tapping the wall of the container 20 times and the height of the column was measured. In this way, we obtained granular beds with reproducible bulk solid volume fractions $\phi \approx 0.54$ (Fig. 1b). After a granular packing was prepared, a fluidization or defluidization experiment was started.

In the fluidization experiments, the air flow was gradually increased and the fluid pressure was monitored. To ensure that a steady-state fluid pressure was established we waited several minutes before each pressure measurement. For each $U$, the generated fluid pressure was measured over time by the pressure sensors that each took 500 samples with a sampling rate of $\approx 41$ Hz. Next, for each sensor, a time averaged fluid pressure $P$ was calculated by averaging over the 500 samples ($\approx 12$ s).

In defluidization experiments, to obtain a well defined initial pressure, the pressure measurements started several seconds before the airflow was stopped. A solenoid valve (V) placed at the entrance of the wind-box was then used to almost instantaneously ($\approx 0.01$ s) stop the air flow. To accurately determine the moment the air flow was stopped, corresponding to $t = 0$ in the experiments, an electric signal generated upon closing the valve was used to trigger the acquisition of fluid pressure. The fluid pressure in time was acquired by taking more than 2000 samples at a sampling rate of about 42 Hz. During these experiments the height of the column was monitored using a ruler, and we estimated change of the height to be less than 0.1%.
FLUIDIZATION AND MODEL

We start by investigating the fluidization of granular columns. For this, reproducible granular beds of height $H$ are prepared by loading a known mass of particles ($d = 80$ µm) into the container from the top (for details, see methods section). Next, an upward vertical air flow is slowly increased until the desired velocity $U$ is reached and then the generated fluid pressure $P$ in the column is measured at three different positions by the differential pressure sensors $S_1$, $S_2$ and $S_3$ (Fig. 1). We test eight different column heights ranging from ≈ 0.4 – 2 m, and for each height we repeat the experiment three times. The pore fluid pressure is then averaged over the three experiments and evaluated as a function of $U$ (Fig. 2a). In general, we observe the well known features of granular column fluidization [1, 2], namely: (i) a first regime where the pore pressure increases linearly with $U$, and, (ii) a second regime where the pore pressure reaches a plateau value $P_m$ independent of $U$. The linear increase of $P_m$ with $H$ is well predicted by the lithostatic pressure $P_\ell(z) = \phi \rho_s g (H - z)$ (see Appendix A). The transition between the two regimes occurs at $U_{mf}$.

In a simple configuration such as that shown in figure 1, Darcy's equation (1) can be used to explicitly deduce the pressure $P$ by noting that the pressure at the top corresponding to $z = H$ is given by $P(H) = 0$ (assuming that the pressure at the top of the column equals the atmospheric pressure). We deduce

$$P(z) = \frac{\mu}{k} U(H - z). \quad (3)$$

This yields that the fluid pressure drop $P(z_i)/(H-z_i)$ obtained at different heights $z_i$, $i \in \{1, 2, 3\}$, is proportional to $U$ and is independent on $H$. However, experimentally, it is clear that the pressure drop decreases with the height of the granular column (Fig. 2 inset). In the following analysis we excluded only the data for $U = 4.7$ mm.s$^{-1}$ for which nonlinear effects in the granular column —such as air bubbles and material dilation— were observed.

Presenting the fluid pressure data as a function of depth in the granular column for different air velocities (Fig. 2d), highlights a nonlinear fluid pressure profile within the granular bed, which indicates non-homogeneous fluidization. This suggests that there exists no unique value of $U_{mf}$ to define the global state of fluidization of the bed (see Appendix A). Similarly, Delebarre [3] argued that, for long columns, $U_{mf}$ cannot be considered as a criterion to characterize powders. To capture the entire pressure profiles observed in our experiments, we propose the following empirical relation (Fig. 2d):

$$P(z) = aU \frac{1 - e^{-b(H-z)}}{b}. \quad (4)$$

Here $a = 4 \times 10^6$ Pa.s.m$^{-2}$ and $b = 0.3$ m$^{-1}$ are global pa-

<figure><img src="image_3.png" /><figcaption>FIG. 2. Pore fluid pressure $P$ versus the air velocity $U$ for different $H$. (a) $P$ measured by sensor $S_1$ at position $z_1 = 0.04$ m. (b) $P$ measured by sensor $S_2$ at position $z_2 = 0.13$ m. (c) $P$ measured by sensor $S_3$ at position $z_3 = 0.3$ m. Insets are the ratio $P/(H - z_i)$ versus $U$, showing a clear decrease of the slope for increasing $H$. Error bars representing the standard error of the mean, obtained over three experiments, are nearly equal or smaller than the symbols. Black stars indicate $U_{mf}$ estimated from the intersect of two linear fits (see Appendix A). (d) Fluid pressure data presented in (a-c) as a function of depth in the granular bed. For $U \in \{0.9, 1.9, 2.4, 2.8, 3.3, 3.8, 4.2\}$ mm.s$^{-1}$, corresponding to symbols circle, diamond, downward triangle, square, plus, triangle and x-mark, respectively. Solid lines are fits by eq. (4). The red dashed line is $P_\ell$ calculated from the bulk $\phi$.</figcaption></figure>

rameters obtained from fitting $P(z)$ across a wide range of $U$ (Fig. 2d). Note that $a$ corresponds to the parameter $\mu/k$ of Darcy's law. For our glass particles, $\mu/k \approx 2.2 \times 10^6$, which agrees reasonably well with $a$ obtained from the fits. The parameter $1/b$ can be interpreted as a typical length scale.

The choice of an approximation given by eq. (4) is no accident and has several advantages:

1. It provides the pressure profile over the entire column and shows that the degree of fluidization depends on the position in the column.

2. For small column height, that is, for small values of $H - z$, the linear nature of the Darcy law, given by the profile (3), is restored.

3. The pressure given by eq. (4) is the solution of a simple differential equation generalizing the usual Darcy law (1):

$$\frac{\partial P}{\partial z} = -aU + bP. \quad (5)$$
4. As we shall see later, this model also helps us to understand the defluidization phenomenon.

Although we cannot identify the origin of the nonlinear pressure profile, we note that the inertial term in eq. (2) should not play a major role. Indeed, experimental conditions indicate that the viscous term in eq. (2) is 10 times greater than the inertial term, consistent with modest particle Reynolds numbers ($Re_p \le 0.02$). For comparison, the pre-Darcy regime occurs at significantly lower $Re_p \approx 10^{-5}$ [13], while the upper limit for the applicability of Darcy's law is typically considered to be $Re_p \sim 1-10$ [14].

Note that eq. (4) is reminiscent of the Janssen model [15], which could thus be a physical phenomenon involved. However, the characteristic length scale $(1/b) \approx 3$ m is nearly 50 times larger than the Janssen length scale $\lambda = R/(2\mu_w K) \sim 0.14$ m, with radius of the container $R$, the effective wall-particle friction coefficient $\mu_w$ and the earth pressure coefficient $K$ ($\mu_w = 0.5$, $K = 0.5$ are typical values for glass bead assemblies [16]). This suggests that $K \mu_w$ is significantly reduced in our experiments. Since $K = 0$ corresponds to a solid and $K = 1$ to a fluid, we expect $K$ to approach 1 in the fluidized beds, implying a significant decrease in $\mu_w$. Although an important reduction in $\mu_w$ has been observed in particle suspensions in water due to lubrication effects [16, 17], this effect seems less likely to occur in air because of its different properties, such as lower viscosity and density.

While we cannot rule out the Janssen effect, we propose that other physical phenomena may also contribute. For example, airflow can alter the particle contact network, potentially through mechanisms such as buckling of force chains [18], leading to hardening of the granular column without changes in packing density. For our glass beads and typical confining pressures (1-10 kPa), Hertz law estimates elastic deformations of the beads to be around 1-5 nm, several orders of magnitude smaller than the grain size. As a result, even tiny deformations can alter the contact network and modify bulk properties, such as force chain lengths and anisotropy, without significantly affecting the macroscopic packing density [19, 20]. Consistent with this idea, a previous study has shown that water flow through a confined assembly of photoelastic disks increases the contact force network [21].

## DEFLUIDIZATION

To test our model, we now investigate the fluid pressure decay over time. For this, a granular column is first fluidized and subsequently, the air supply is instantly stopped and the decay of fluid pressure is measured in time. Similar to the fluidization experiments discussed before, we vary the height of the column and for each

height the experiment is repeated three times. During these experiments, we do not observe any visible rearrangement of the grains nor do we measure any significant change in the height of the granular columns. In general, we find that once the air supply is stopped the pore fluid pressure in the non-expanded columns decreases in time, first rapidly and then at a slower rate (Fig. 3). The initial fluid pressure at $t = 0$ as a function of position in the granular bed is well predicted by eq. (4) (Fig. 4a) with $a$ and $b$ deduced from the fluidization experiments.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Representative experiments showing the pore fluid pressure decrease in time, measured using the three pressure sensors for granular columns of increasing height $H$ (from dark purple to light pink and indicated by the arrow). (a) For $U = 1.4$ mm.s$^{-1}$. (b) $U = 2.4$ mm.s$^{-1}$. (c) $U = 2.8$ mm.s$^{-1}$. (d) $U = 3.3$ mm.s$^{-1}$. Solid black lines are fits of eq. (9).</figcaption></figure>

The relation for the time-dependent pressure evolution is obtained from the local mass conservation law for the fluid in the column. This relation involves the apparent density $(1 - \phi)\rho$ of the gas, and the microscopic velocity of the gas $u$. It writes:

$$ \frac{\partial}{\partial t}((1 - \phi)\rho) + \frac{\partial}{\partial z}((1 - \phi)\rho u) = 0. \quad (6) $$

The link between the microscopic velocity $u$ and the superficial fluid velocity $U$ being given by $U = (1 - \phi)u$, we deduce

$$ \frac{\partial}{\partial t}((1 - \phi)\rho) + \frac{\partial}{\partial z}(\rho U) = 0. \quad (7) $$

Assuming an ideal gas, that is the total pressure is proportional to the density: $P_{atm} + P = \beta\rho$ and the experimental constant $\beta = RT/M$ (Pa.m$^3$.kg$^{-1}$), with the molecular mass of air $M$, gas constant $R$ and absolute temperature $T$. We consider that the solid volume fraction $\phi$ is constant and using the extended Darcy's law (5)
derive the law for the evolution of the fluid pressure:

$$ \frac{\partial P}{\partial t} + \frac{1}{(1-\phi)a} \frac{\partial}{\partial z} \left( (P_{\text{atm}} + P)(bP - \frac{\partial P}{\partial z}) \right) = 0. \quad (8) $$

It is therefore relevant to assume that the over-pressure P satisfies $P \ll P_{atm}$ and neglect some terms (see Appendix B for further details and dimensional analysis). We obtain

$$ \frac{\partial P}{\partial t} - \frac{P_{\text{atm}}}{(1-\phi)a} \frac{\partial^2 P}{\partial z^2} + \frac{P_{\text{atm}}b}{(1-\phi)a} \frac{\partial P}{\partial z} = 0. \quad (9) $$

Compared with the usual diffusion model based on Darcy's law (1), which formally corresponds to the case where $b = 0$, we obtain an equation containing an additive first-order term. The experimental pressure evolution is fitted by numerically solving eq. (9) using MATLAB's lsqcurvefit and ODE45 functions, and $\phi$ is the only free fitting parameter. The initial condition is given by eq. (4) corresponding to the steady-state solution for a given velocity $U$, that is

$$ P(0, z) = aU \frac{1 - e^{-b(H-z)}}{b}. \quad (10) $$

By imposing the following conditions at the top ($z = H$) and bottom ($z = 0$) of the domain

$$ P(t, H) = 0 \quad \text{and} \quad \frac{\partial P}{\partial z}(t, 0) = bP(t, 0), \quad (11) $$

we simulate the stopping of the air injection starting from $t = 0$ (since $U = 0$ corresponds to $\partial_z P = bP$ in the generalized Darcy law (5)).

We find good agreement between the fits and the experimental defluidization over time for different $z$, $H$ and $U$ (Fig. 3). The inferred $\phi$, shown in Fig. 4b, is generally close to the experimentally measured bulk value $\phi \approx 0.54$, and only for the smallest $H-z$ do we observe deviations. This is likely related to $a$ and $b$, which were obtained from global fits that include fluidization experiments at higher $U$, where nonlinear phenomena, like particle sliding at the contact scale or the emergence of tiny air bubbles, may come into play.

Importantly, we note that, unlike the diffusion model (see Appendix C), our extended Darcy law-based approach yields a $\phi$ that appears independent on $U$ and $H-z$.

## CONCLUSION

Our findings reveal that the fluid pressure profile in fluidized granular columns is nonlinear. This behavior can be captured by an empirical term in Darcy's law. The key advantage of this extended Darcy law is that

<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Initial fluid pressure at $t = 0$ s as a function of position in the granular bed averaged over three defluidization experiments for $U = 1.4, 2.4, 2.8, 3.3 \text{ mm.s}^{-1}$ (circle, downward triangle, square, triangle). Error bars representing the standard deviation over three experiments are nearly equal or smaller than the symbols. Black dashed lines are eq. (5) with $a$ and $b$ from the fluidization experiments. (b) Solid volume fraction $\phi$ obtained from fitting the defluidization $P(t)$ by eq. (9), as a function of position in the column. The dashed line represents $\phi = 0.54$.</figcaption></figure>

it can be easily rewritten and used in the derivation of the pressure evolution law unifying the fluidization and defluidization processes in granular columns.

<footnote>* siet.van_den_wildenberg@uca.fr</footnote>

<footnote>[1] O. Roche, Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective, Bull Volcanol **74**, 1807 (2012).</footnote>

<footnote>[2] D. Geldart, *Gas fluidization Technology* (Wiley NY, 1986).</footnote>

<footnote>[3] A. Delebarre, Does the minimum fluidization exist?, *Journal of Fluids Engineering* **124**, 595 (2002).</footnote>

<footnote>[4] S. Ergun, Fluid flow through packed columns, *Chemical Engineering Progress* **48**, 89 (1952).</footnote>

<footnote>[5] A. Anantharaman, R. A. Cocco, and J. W. Chew, Evaluation of correlations for the minimum fluidization velocity (umf) in gas-solid fluidization, *Powder Technology* **323**, 454 (2018).</footnote>

<footnote>[6] M. Rasteh, G. Ahmadi, and S. H. Hosseini, Effect of gaus- sian size distribution owidth on minimum fluidization ve- locity in tapered gas-solid fluidized beds, *Particuology* **66**, 71 (2022).</footnote>

<footnote>[7] A. Rao, J. S. Curtis, B. C. Hancock, and W. C., The effect of column diameter and bed height on minimum fluidization velocity, *American Institure of Chemical En- gineers* **56**, 2304 (2010).</footnote>

<footnote>[8] H. Quan, N. Fatah, and C. Hu, Diagnosis of hydrody- namic regimes from large to micofluidized beds, *Chemical Engineering Journal* **391**, 123615 (2020).</footnote>

<footnote>[9] Z. Jiang and N. Fatah, New investigation of micro- fluidized bed: The effect of wall roughness and parti- cle size on hydrodynamics regimes, *Chemical Engineering Journal* **430**, 133075 (2022).</footnote>

<footnote>[10] E. Breard, L. Jones, J. and Fullard, G. Lube, C. Davies, and J. Dufek, The permeability of volcanic mixtures- implications for pyroclastic currents, *JGR solid earth* **10**, 1343 (2018).</footnote>

<footnote>[11] S. Leibrandt and J. L. Le Pennec, Towards fast and routine analyses of volcanic ash morphometry for erup-</footnote>
tion surveillance applications, Journal of Volcanology and Geothermal Research 297, 11 (2015).

[12] E. C. P. Breard, J. Dufek, and G. Lube, Enhanced mobility in concentrated pyroclastic density currents: An examination of a self-fluidization mechanism, Geophysical Research Letters 45, 654 (2018).

[13] R. M. Fand, B. Y. K. Kim, A. C. C. Lam, and R. T. Phan, Resistance to the flow of fluids through simple and complex porous media whose matrices are composed of randomly packed spheres, J. Fluids Eng. 109, 268 (1987).

[14] W. Sobieski and A. Trykozko, Darcy's and forchheimer's laws in practice. part 1. the experiment, Technical Sciences 17, 321 (2014).

[15] H. A. Janssen, Versuche über getreidetdruck in silozellen, Z. Ver. Dtsh. Ing 39, 1045 (1895).

[16] S. M. Brice, M. Georgelin, S. Deville, and A. Pocheau, Wall friction and janssen effect in the solidification of suspensions, Soft Matter 14, 9498 (2018).

[17] I. Rosenhek-Goldian, N. Kampf, and J. Klein, Trapped aqueous films lubricate highly hydrophobic surfaces, ACS Nano 12, 10075-10083 (2018).

[18] J. D. Goddard, Nonlinear elasticity and pressure-dependent wave speeds in granular media, Proc. R. Soc. London, Ser. A 430, 105 (1990).

[19] C. F. Schreck, T. Bertrand, O. C. S., and M. D. Shattuck, Repulsive contact interactions make jammed particulate systems inherently nonharmonic, Phys. Rev. Lett. 107, 078301 (2011).

[20] E. T. Owens and K. E. DANIELS, Sound propagation and force chains in granular materials, EPL 94, 54005 (2011).

[21] N. Mahabadi and J. Jang, The impact of fluid flow on force chains in granular media, Appl. Phys. Lett. 110, 041907 (2017).

<figure><img src="image_6.png" /><figcaption>FIG. 5. Pore fluid pressure $P$ versus the air velocity $U$ for different column heights $H$, same data as Fig. 2 in the main text.</figcaption></figure>

(a,b and c) $P$ measured by sensor $S_1$, $S_2$ and $S_3$, respectively, as a function of $U$. Dashed lines are linear fits to the data in the regime where the pore pressure increases with $U$. Dotted lines are estimations of the plateau pressure. (d) Maximum fluid pressure $P_m$ obtained from the plateaus in (a-c), versus $H-z$. Colors correspond to the colors (a-c), and the black line indicates the lithostatic pressure, with the bulk solid volume fraction $\phi=0.54$. (e) Minimum fluidization velocity $U_{mf}$ estimated from the intersection between the dashed and dotted lines in (a-c) as a function of $H-z$.

## B - DIMENSIONAL STUDY

In this appendix, we explain the transition from equation (8) to its simplified version given by (9). The principle is to justify that certain terms in equation (8) are negligible compared to other terms. In order to compare the terms, we must first write a dimensionless version of equation (8). To do this, we introduce the following dimensionless variables and unknowns - which will be noted with a symbol $\sim$:

$$t = T\tilde{t}, \quad z = H\tilde{z}, \quad P = \phi\rho_s g H \tilde{P}.$$

where $H$ stands for the characteristic length of the experiment, and $T$ for a characteristic time. The equation (8) can therefore be written as

$$\frac{\partial \tilde{P}}{\partial \tilde{t}} = \tilde{D} \frac{\partial}{\partial \tilde{z}} \left( \frac{\partial \tilde{P}}{\partial \tilde{z}} - \tilde{b} \tilde{P} + \tilde{\varepsilon} \tilde{P} \frac{\partial \tilde{P}}{\partial \tilde{z}} - \tilde{\varepsilon} \tilde{b} \tilde{P}^2 \right), \quad (12)$$

where $\tilde{D}$, $\tilde{b}$ and $\tilde{\varepsilon}$ are three dimensionless numbers given by

$$\tilde{D} = \frac{T p_{\text{atm}}}{(1 - \phi) a H^2}, \quad \tilde{b} = bH \quad \text{and} \quad \tilde{\varepsilon} = \frac{\phi \rho_s g H}{p_{\text{atm}}}.$$

Thus, depending on the orders of magnitude of the coefficients $b$ and $\tilde{\varepsilon}$, several regimes can be distinguished. Knowing the values $\phi = 0.54$, $\rho_s \approx 2500 \text{ kg.m}^{-3}$, $p_{\text{atm}} = 10^5 \text{ Pa}$, $b = 0.3 \text{ m}^{-1}$ and $g = 9.8 \text{ m.s}^{-1}$, we can separate the cases simply as a function of the height $H$ of the column:

# APPENDIX

## A - ESTIMATION $U_{mf}$

For each column height, the minimum fluidization velocity is determined by identifying the intersection point between a linear fit of the data in the regime, where $P$ increases with $U$, and a fit of the pore pressure in the plateau region. The inferred $U_{mf}$ increases with $H$ as observed in [7].
1. Case $H = 5$ cm:

$\tilde{b} \approx 1.5 \times 10^{-2} \ll 1$ and $\tilde{\epsilon} \approx 6.6 \times 10^{-3} \ll 1$ so that (12) can be reasonably approximated by a simple diffusion equation:

$$ \frac{\partial \tilde{P}}{\partial \tilde{t}} = \tilde{D} \frac{\partial^2 \tilde{P}}{\partial \tilde{z}^2}. \qquad (13) $$

In this case, the initial and boundary conditions (10) and (11) can be approximated by

$$ \tilde{P}(0, \tilde{z}) = \frac{aU}{\phi \rho_s g} (1 - \tilde{z}) \qquad (14) $$
and $\tilde{P}(\tilde{t}, 1) = \frac{\partial \tilde{P}}{\partial \tilde{z}}(\tilde{t}, 0) = 0$

so that an explicit solution can be easily obtained by separating the variables and using Fourier series:

$$ \tilde{P}(\tilde{t}, \tilde{z}) = \sum_{n=0}^{+\infty} a_n \cos\left(\frac{(2n+1)\pi}{2}\tilde{z}\right) \exp\left(-\left(\frac{(2n+1)\pi}{2}\right)^2 \tilde{D}\tilde{t}\right), \qquad (15) $$

where the coefficients $a_n$, $n \ge 0$, are given using the initial data:

$$ a_n = \frac{8aU}{\phi \rho_s g (2n+1)^2 \pi^2}. \qquad (16) $$

2. Case $H = 40$ cm (order of magnitude of experiments carried out):

$\tilde{b} \approx 1.2 \times 10^{-1}$ and $\tilde{\epsilon} \approx 5.3 \times 10^{-2} \ll 1$ so that (12) can be approximated by the equation (9) (written here in dimensionless form):

$$ \frac{\partial \tilde{P}}{\partial \tilde{t}} = \tilde{D} \frac{\partial^2 \tilde{P}}{\partial \tilde{z}^2} - \tilde{D} \tilde{b} \frac{\partial \tilde{P}}{\partial \tilde{z}}. \qquad (17) $$

Note that, although equation (17) remains linear, the boundary conditions (11) do not allow us to obtain a simple analytic expression (using Fourier series, for example) for the solution.

3. Case $H = 5$ m:

$\tilde{b} \approx 1.5$ and $\tilde{\epsilon} \approx 6.6 \times 10^{-1}$ and no term in (8) seems negligible.

## C - COMPARISON WITH DARCY'S MODEL

In this section, we compare our model in eq. (9) for the defluidization of a granular column to a model based on Darcy's law eq. (1). Following the same procedures as in the main text and Appendix B, we derive the 1D-diffusion equation [1]. In order to complete the model, we need to impose initial conditions corresponding to the situation at the beginning of the experiment, and boundary conditions indicating interactions with the outside. We

start with a partially fluidized condition (a proportion $\zeta \in [0, 1]$ of the complete fluidization condition) and assume that at the top of the column the pressure equals the atmospheric pressure, while at the bottom of the column a zero air injection velocity $U$ is imposed:

$$ \left\{ \begin{array}{ll} \frac{\partial P}{\partial t} = \mathcal{D} \frac{\partial^2 P}{\partial z^2} & \text{for } t > 0 \text{ and } z \in (0, H) \\ P(0, z) = \zeta \phi \rho_s g (H - z) & \text{for } z \in (0, H) \\ P(t, H) = \frac{\partial P}{\partial z}(t, 0) = 0 & \text{for } t > 0 \end{array} \right. \qquad (18) $$

The pressure evolution is fitted by numerically solving the coupled eq. (18) using MATLAB's lsqcurvefit and ODE45 functions, with two free fitting parameter $\zeta$ and $\mathcal{D}$. The fits to the data as well as the inferred parameters are shown in Fig. 6. Generally, we find that both models fit the data similarly well, although the diffusion model uses two parameter optimization in contrast to the single parameter in eq. (9). The inferred $\zeta$ decreases with $H-z$ and increases with $U$. Importantly, the inferred $\mathcal{D}$ increases about 40% for the $H-z$ tested, consistent with previous observations [1, 12].

<figure><img src="image_7.png" /><figcaption>FIG. 6. Fluid pressure versus time for different column heights, same experimental data as in Fig.3. (a) for $U = 2.4 \text{ mm.s}^{-1}$. (b) $U = 3.3 \text{ mm.s}^{-1}$. Black lines represent fits of (18). For comparison, the fits using the extended Darcy law model are also shown (gray dashed lines). Both models describe the data for all H and U values equally well. (c) Degree of fluidization $\zeta$ obtained from the fits decreases with increasing column height $H$ but increases with $U$, for $U \in \{2.4, 3.3\} \text{ mm.s}^{-1}$, corresponding to symbols downward triangle and triangle, respectively. (d) Diffusion constant $\mathcal{D}$, obtained from the fits, shown as a function of $H$. For comparison, values of $\phi$ from Fig.4b, are shown in light gray.</figcaption></figure>